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Abstract 


Stress  estimates  obtained  from  modelling  various 
regions  near  the  Middle  America  trench  have  been  combined 
with  rock  failure  criteria  expressed  in  statistical  terms  to 
yield  probabilistic  measures  of  seismic  risk  in  these  areas. 
These  measures  allow  the  study  of  different  geodynamic 
processes,  in  such  a  way  that  the  probable  location  of 
seismic  failure  can  be  estimated.  Repetitive  modelling  of  a 
system  evolving  in  time  can  suggest  the  most  likely  time  of 
failure . 


A  new  measure  of  seismic  instability  has  been  developed 
which  combines  the  stress  components  into  a  unique  measure 
of  seismic  risk.  This  function  depends  on  the  material 
properties,  and  the  geological  structures  involved  in  a 
particular  problem. 

The  results  obtained  here  for  the  Valley  of  Mexico  are 
in  agreement  with  the  observed  seismicity  and  the  geodynamic 
knowledge  of  the  area.  This  seismic  activity  seems  to  be 
triggered  by  the  reduction  of  water  in  the  aquifer  that 
underlies  Mexico  City. 

The  distance  at  which  the  stress  drop  due  to  the 
occurrence  of  an  earthquake  can  affect  the  seismic  behaviour 
of  distant  locations  is  linked  to  the  thickness  of  the 
elastic  part  of  the  tectonic  plate  at  the  epicentral 
location.  This  variation  on  the  possible  range  of  seismic 
activity  supports  the  existence  of  long  term  precursors  for 
distances  up  to  103km. 
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The  high  risk  areas  derived  from  the  models  can  be 
identified  with  seismic  active  areas.  This  is  done  with 
simple  models  of  the  geodynamic  processes  involved.  The 
material  properties  used  are  in  general  the  elastic 
parameters  and  their  variation  with  temperature.  The 
external  forces  applied  to  the  models  are  those  that  have 
the  potential  to  change  the  overall  distribution  of  seismic 
instability  such  as  the  gravitational  pull  of  the  downgoing 
slab,  the  viscous  drag  of  the  mantle  and  the  frictional 
force  between  plates.  Therefore,  these  distributions  are 
regional  in  nature. 

In  spite  of  these  simplifications,  the  distributions  of 
seismic  instability  calculated  for  the  Middle  America  trench 
region  are  in  accordance  with  the  observed  seismic  behaviour 
of  the  area.  Shallow  activity  (depth<70  km)  along  both  sides 
of  the  trench  is  more  frequent  than  in  the  North 
Amer ica-Car ibbean  plate  boundary.  Deep  earthquakes  are  more 
likely  to  occur  in  the  trench  south  of  the  Tehuantepec 
ridge.  Events  in  this  southern  branch  have  a  smaller  range 
of  effects  than  those  to  the  north  of  the  ridge  due  to  the 
difference  in  the  cohesion  between  the  plates  and  the 
subducting  angle. 
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1  .  Introduction . 


'Earthquakes  can  be  treated  as  mechanical  instabilities,, 
They  occur  when  the  rock  is  unable  to  support  an  increase  of 
stress.  Therefore,  seismic  risk  studies  are  approached  here 
by  determining  how  stable  the  rocks  are  under  various  stress 
conditions . 

Many  failure  criteria  exist  for  rocks  and  it  is  not 
obvious  that  any  of  them  are  representative  for  the 
geological  scales  used  here.  Nevertheless,  it  is  possible  to 
overcome  this  problem  by  using  a  generalized  definition  of 
seismic  instability  independent  of  the  type  of  failure 
envelope  used.  This  will  allow  the  failure  criteria  to  be 
changed  without  losing  any  generality. 

Different  geodynamic  processes  will  be  treated 
throughout  this  work  in  which  it  will  be  necessary  to  vary 
the  material  properties  and  the  stress-strain  relations  that 
describe  the  behaviour  of  the  materials  involved.  These 
different  exercises  will  provide  an  insight  into  the 
geodynamical  processes  occurring  in  the  Middle  America 
trench  region,  as  well  as  illustrate  what  kind  of  studies 
can  be  done  with  these  kinds  of  procedures. 

The  results  obtained  can  be  generalized  to  other  areas. 
Explanations  of  some  observed  seismic  phenomena  can  be 
derived,  (eg.  Why  does  the  range  of  long  term  precursors 
vary?).  Instability  distributions  will  be  shown  as 
predictors  of  the  possible  location  of  failure  initiation  or 
as  time  predicting  functions  for  particular  locations. 
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Individual  studies  of  different  seismically  active 
areas  related  geodynamically  to  the  different  parts  of  the 
Middle  America  trench  give  an  overall  picture  of  the  region. 
Stresses  calculated  by  means  of  the  finite  elements 
technique  will  be  used  in  some  cases,  analytical  methods 
will  be  applied  in  others. 

The  determination  of  the  stability  of  state  of  stress 
in  the  earth  is  of  vital  importance  in  understanding  any  of 
the  different  processes  which  occur  in  the  interior  of  the 
earth.  In  this  chapter  the  current  knowledge  of  the 
propagation  of  stresses  in  the  upper  part  of  the  earth  is 
reviewed.  I  briefly  mention  the  different  techniques  that 
exist  for  measuring  these  stresses,  I  explore  the  possible 
relation  between  earthquakes  and  stresses;  between  the 
tectonic  models  of  the  earth  and  the  location  of  the 
observed  seismic  activity;  and  the  consequences  that  the 
different  kinds  of  seismic  behaviour  observed  might  have  in 
our  understanding  of  the  earth. 

In  chapter  2,  I  connect  these  stress  arguments  to  the 
occurrence  of  seismic  failures  by  means  of  a  definition  of 
seismic  stability.  Once  defined  this  concept  can  undergo 
several  variations  which  are  discussed  in  sections  2.3  to 
2.6.  These  individual  forms  are  applied  throughout  chapters 
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1.1  The  Problem. 

Geodynamics  is  the  quantitative  study  of  the  processes 
that  take  place  in  the  real  earth.  A  major  problem  that  it 
has  to  overcome  is  that  direct  observations  of  these 
processes  in  the  interior  of  the  earth  are  impossible; 
knowledge  relies  on  the  interpretation  of  a  series  of 
geophysical  observables,  from  which  reasonable  deductions 
may  be  made.  This  thesis  deals  with  the  geodynamical 
mechanical  state  of  the  earth  and  this  mechanical  state  can 
be  studied  through  the  distribution  and  propagation  of 
stresses  and  strains  throughout  the  different  materials  that 
constitute  the  earth.  In  particular,  the  way  the  stress 
state  may  be  related  to  seismic  instability  is  examined 
here . 

1.2  Stress  and  Strain  in  the  Earth  Sciences. 

The  concepts  of  stress  and  strain  are  fundamental  to 
the  interpretation  of  earthquake  mechanisms,  and  any  study 
of  structural  geology  needs  them  for  any  description  of  rock 
behaviour.  The  stresses  that  exist  in  a  given  rock  mass  or 
geological  area  control  the  processes  that  result  in  the 
mechanical  state  of  this  material,  and  allow  us  to  make 
predictions  of  the  changes  that  might  occur  in  the  future. 

Strains  and  displacements  are  frequently  measured,  and 
from  these  the  stresses  can  be  deduced.  The  stresses  at  a 
point  in  the  interior  of  a  body  are  determined  by  the  system 
of  forces  that  act  in  the  vicinity  of  that  given  point.  The 
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deformation  of  the  body  in  the  vicinity  of  such  point  is 
termed  strain  (Richter,  1958).  The  analysis  of  strain  is 
fundamental  for  studies  of  the  deformation  of  any  material 
(Jaeger  and  Cook,  1979).  The  study  of  stress  is  in  itself  a 
pure  static  analysis,  independent  of  the  properties  assumed 
for  the  material  which  may  be  elastic,  plastic  or  any  other. 
This  means  that  in  order  to  state  the  amount  of  stress  on  a 
given  point  inside  a  body  it  is  not  necessary  to  specify  the 
physical  properties  of  the  body  itself. 

Knowledge  of  the  history  of  the  stress  state  is  also 
fundamental  to  the  understanding  of  the  tectonic  processes 
that  occur  in  the  earth.  However  stress  values  are  very 
difficult  to  obtain  even  for  a  single  location  of  the  earth 
lithosphere;  Wyss  (1977)  mentions  this  difficulty  and  points 
out  that  of  the  six  components  of  the  stress  tensor  usually 
only  a  few  can  be  estimated. 

Rock  behaviour  under  stress  changes  is  poorly 
understood,  since  it  depends  on  many  factors  like  the  forces 
that  have  been  applied  previously,  the  duration  of  such 
forces,  the  confining  pressure  and  the  temperature.  In 
geodynamical  terms  the  state  of  stress  is  linked  to  causes 
(loading  and  unloading,  heating  and  cooling,  plate  motions 
and  driving  forces,  etc.),  to  consequences  (creep 
deformation,  and  seismic  failure),  and  to  rheology  (the 
depth  over  which  stress  can  be  supported  and  the  time 
dependence  of  material  properties).  None  of  these  links  has 
been  characterized  in  sufficient  detail  to  define  the 
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constitutive  equation  for  any  part  of  the  lithosphere 
without  making  many  assumptions  (Solomon  et  al.,  1980). 

Due  to  the  lack  of  understanding  of  the  earth  processes 
and  the  relation  between  the  strains  and  stresses  at  depth, 
the  constitutive  relationship  is  poorly  understood.  Several 
ways  to  relate  them  have  been  proposed;  it  has  been  found 
from  laboratory  tests  that  strain  in  rocks  increases  in  a 
linear  way  when  they  are  subjected  to  small  external 
stresses.  This  kind  of  linear  dependence  between  stress  and 
strain  is  named  linear  elasticity  and  is  characterized  by 
the  generalized  Hooke’s  law 


where  o  j  j  are  the 
C  j  j  m  n  the  elastic 
the  infinitesimal 


cartesian  components  of  the  stress  tensor, 
modulus  components.  The  components  emn  of 
strain  tensor  are  defined  in  terms  of  u  by 


£  m  n  ;  m^U|],(  n  )/2 


(1.2) 


where  u(x)  is  the  distortion  suffered  by  a  body  at  point  x 
when  a  set  of  forces  is  applied  to  it,  the  displacement  of  a 
particle  originally  at  x,  from  the  initial  state  to  the 
distorted  one.  The  initial  state  need  not  be  stress  free. 

While  this  kind  of  constitutive  relation  seems  adequate 
for  many  geophysical  and  engineering  purposes,  there  is  need 
for  a  more  sophisticated  formulation  to  account  for  the 
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non-linear  behaviour  observed  in  rocks  when  the  temperature 
and  confining  pressure  increase.  One  equation  that  seems  to 
describe  these  observations  is  that  of  a  viscous  or 
visco-elastic  material  (Jaeger  and  Cook,  1979), 

o  © 

e  =E '  1  a  +o/r}  (1.3) 

where  E  is  the  Young's  modulus  and  77  the  viscosity  of  the 
material. 

The  dependence  of  equation  1.3  on  temperature  and 
confining  pressure  is  implied  in  the  viscosity.  It  is 
important  to  note  that  this  equation  predicts  that: 

1.  For  high  stress  rate  ( o  )  the  first  term  will  dominate 
and  the  material  will  behave  elastically. 

2.  For  small  o  values  the  second  term  will  dominate  giving 
the  material  a  viscous  behaviour. 

3.  As  the  confining  pressure  and  the  temperature  increase, 
the  material  becomes  gradually  less  viscous  until  it 
reaches  a  "fluid"  state. 

Savage  and  Prescott  (1978),  proposed  this  kind  of 
behaviour  to  model  depth  dependence  of  the  stresses  in  the 
earth.  Ringwood  (1975),  also  suggested  a  non-linear  kind  of 
stress  variation  with  depth  for  ocean  and  continental  rocks. 
Kirby  (1977),  explains  that  in  experiments  conducted  in 
cylinders  of  rock  subject  to  a  uniaxial  stress  state 
a i i 2 2=o 3 3=P  where  the  x,  direction  is  parallel  to  the 
cylinder  axis  and  P  is  the  confining  pressure,  the 
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differential  stress  a=|a11-P|  is  related  to  the  strain  rate 
e  ,  temperature  T  and  pressure  P  for  most  earth  materials  by 
the  non-linear  expression 

e  =A ' [ s inh ( aa ] n exp[ - ( E ' +PV' ) /RT]  (1.4) 

where  A’,E',a,7j  and  V'  are  material  properties  and  R  is  the 
gas  constant.  Equation  1.4  for  large  temperature  and 
pressure  becomes  the  exponential  law 

e  =A  exp(/3  o)  exp[  -  ( E ?  +PV’  )  /RT  ]  (1.5) 

where  A=A'/2  and  /3  =  77a. .  This  type  of  equation  (1.5)  is  taken 
by  Kirby  (1977)  and  Savage  and  Prescott  (1978)  and  others, 
as  representative  of  the  upper  part  of  the  mantle.  In  this 
thesis  the  elastic  behaviour  of  the  upper  part  of  the  crust 
and  the  viscoelastic  nonlinear  behaviour  of  the  stresses  in 
depth  and  temperature  will  be  used  extensively. 

It  seems  probable  that  earthquakes  are  mechanical 
instabilities  evidenced  by  a  sudden  failure  of  the  rock  to 
sustain  the  shear  stress  acting  across  a  surface;  the 
surface  is  usually  an  already  existing  fault  but  may  be  one 
newly  created  by  the  failure.  Due  to  the  large  time  scales 
in  which  geologic  processes  take  place  it  is  possible  to 
picture  the  earth  in  equilibrium  at  any  time.  This  implies 
that  the  rate  of  accumulation  of  stress  in  the  earth’s  upper 
layers  is  relatively  small. 
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Nevertheless,  when  the  accumulated  stress  overcomes  the 
strength  of  the  rocks,  failure  occurs.  Stability  values  can 
measure  the  probability  that  failure  will  occur  in  a 
particular  location  within  a  given  volume.  This  volume  could 
contain  several  materials  and  may  contain  discontinuities. 
The  only  constraint  is  that  the  stress-strain  relationships 
throughout  the  whole  body  have  to  be  known  or  at  least 
assumed . 

Stability  is  then  a  function  of  the  state  of  stress  at 
particular  locations,  it  compares  this  state  of  stress  with 
the  overall  stress,  to  obtain  a  measure  of  how  stable  these 
locations  are  with  respect  to  the  surrounding  areas.  The 
stability  functions  that  have  been  developed  will  be 
discussed  in  chapter  2  of  this  thesis. 

The  reason  why  earthquakes  occur  only  at  particular 
locations  of  the  earth's  lithosphere  is  due  to 
concentrations  of  stresses  at  weakened  areas;  The 

combination  results  in  unstable  areas.  These  concentrations 

* 

are  likely  to  occur  if  forces  are  applied  at  particular 
locations  of  a  body  within  the  neighbourhood  of  the 
application  point.  They  also  appear  in  the  vicinity  of  any 
kind  of  discontinuity,  which  can  be  a  sharp  boundary  or 
other  change  in  physical  properties.  A  simple  example  can  be 
stated.  If  a  force  is  applied  at  the  surface  of  an  elastic 
plate,  the  plate  will  deform  in  the  vicinity  of  the  point  of 
application.  At  locations  away  from  the  force  the  plate  has 
not  suffered  any  distortion  at  all.  This  distortion  is 
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linked  to  the  stresses  by  the  relations  in  section  1.1, 
therefore,  the  stresses  due  to  this  force  diminish  with  the 
distance  to  the  applied  force.  The  rate  of  attenuation  is 
determined  by  the  stress-strain  relationship  characteristic 
for  the  material. 

The  degree  of  instability  is  given  by  the  system  of 
forces  acting  in  the  area  and  the  geological  structures  and 
type  of  material  present  in  it,  and  its  "strength".  All 
these  things  combine  in  such  a  way  that  the  smaller  the 
stability  of  a  location  is,  the  greater  probabilty  it  has  to 
be  the  origin  of  a  seismic  event. 

It  is  the  purpose  of  this  thesis  to  demonstrate  how 
simple  notions  of  rock  mechanics  can  combine  with  what  is 
known  about  earthquakes  to  develop  a  simple  stability 
criterion.  This  criterion  may  give  us  insight  into  aspects 
of  earthquake  generation,  and  probable  causes  for  some  of 
the  earthquake  related  phenomena  observed.  Application  to 
particular  problems  will  be  demonstrated  throughout  the 
following  chapters.  In  order  to  attack  each  of  these 
problems  the  stability  criterion  will  suffer  small  changes 
that  will  be  explained  in  the  corresponding  parts  of  this 
thesis.  However,  the  general  features  of  stability  will 
remain  valid.  In  chapter  2  these  generalities  are  presented. 
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1.3  Stresses  and  Earthquakes. 

An  earthquake  can  be  thought  of  as  a  sudden  release  of 
energy  that  takes  place  at  a  particular  location,  and  is 
associated  with  faulting  or  other  dynamic  processes. 
Hypocenter  is  a  term  used  since  the  1850's  when  Mallet  and 
others  believed  that  an  earthquake  originated  in  a  small 
volume  that  could  be  represented  by  a  point  for  most 
purposes  (Richter,  1958),  and  used  to  designate  the  location 
of  the  origin  of  the  earthquake.  The  term  is  now  thought  of 
as  the  location  of  the  initial  rupture  of  the  rocks,  where 
the  earthquake  originates. 

Richter  (1958),  also  recognized  that  the  energy  source 
for  a  tectonic  earthquake  is  energy  stored  within  the  earth 
during  a  long  period  of  growth  of  strain.  If  the  associated 
stresses  accumulate  continuously  they  will  reach  the  peak 
strength  of  the  rock  and  failure  will  generally  occur.  He 
also  concluded  that  fracture  will  occur  along  preexisting 
weaknesses  and  that  the  active  faults  are  weakness  planes 
where  fractures  occur  one  after  the  other. 

The  seismic  belts  on  the  surface  of  the  earth  are  where 
most  earthquakes  are  located  (figure  1),  (Richter , 1 958 ) . 

This  observation  was  made  by  Gutenberg  (1949),  and  he 
suggested  that  the  stress  in  these  belts  somehow  is 
accumulating  and  as  result  of  this,  earthquakes  are  being 
produced. 

There  are  two  main  types  of  activity  occurring  within 
the  seismic  plate.  One  type  results  from  stresses  due  to  the 
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Figure  1....  seismicity  of  the  earth-seismic  belts. 
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deformation  of  the  plate  itself;  these  are  earthquakes  that 
take  place  close  to  the  plate  boundaries.  They  are  due  to 
the  bending  of  the  plate  (Isacks  et  al.,  1968).  There  are 
also  earthquakes  that  result  from  stresses  transmitted 
through  the  plate  away  from  its  boundaries  (Isacks  and 
Molnar  ,  1971). 

The  implications  of  seismic  observations  are  that  the 
generation  of  earthquakes  in  a  subduction  zone  is  governed 
by  the  distribution  of  stresses  within  the  lithosphere  and 
by  changes  in  the  physical  properties  of  the  material  due  to 
the  presence  of  the  subducting  slab  (Isacks  and 
Molnar , 1 97 1 ) .  Several  hypotheses  exist  to  explain  shear 
fractures  observed  at  intermediate  depths.  The  one  that 
suggests  a  high  pressure,  high  temperature  mechanical 
instability  (Byerlee  and  Brace,  1969)  seems  to  be  the  most 
appropriate.  It  is  expected  that  the  cold  slab  remains 
cooler  than  the  surrounding  medium  (which  was  already  at 
higher  temperatures,  when  the  process  started),  down  to  a 
depth  where  the  slab  cannot  exist  as  a  solid.  This  model  of 
a  cold  subducting  slab  has  been  generaly  accepted  and  the 
temperatures  in  it  have  been  calculated  with  several  models 
(i.e.  McKenzie,  1969;  Minear  and  Toksoz,  1970;  Griggs, 

1972) . 

These  thermal  induced  stresses  have  no  significant 
effect  on  the  mechanisms  of  shallow  earthquakes.  The 
presence  of  this  kind  of  seismic  activity,  may  not  be 
necessarily  related  to  the  plate  boundaries  but  to  local 
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geological  phenomena.  For  the  study  of  this  kind  of 
earthquakes  viscous  behaviour  of  the  surrounding  medium  is 
not  to  be  expected  due  to  the  shallow  depths  at  which  they 
originate. 

1.4  Stress  modelling  in  geodynamics. 

'The  pi ate-tecton ics  hypothesis  explains  the  tectonic 
and  seismic  activity  now  occurring  within  the  upper  layers 
of  the  earth  as  resulting  from  the  interaction  of  a  small 
number  of  large  rigid  plates  whose  boundaries  are  the 
seismic  belts  of  the  world.  The  plates  outlined  by  these 
belts  are  not  actively  deformed  except  along  their  borders 
and  the  motion  occurring  within  them  is  mostly  limited  to 
broad  epirogenic  movements.  These  plates  may  contain 
continental  as  well  as  oceanic  surfaces.  The  seismic  belts 
are  zones  where  differential  movements  between  plates  take 
place'  (Le  Pichon,  1973). 

Plate  boundaries  are  zones  of  concentration  of  stress 
because  they  represent  the  discontinuities  in  the 
lithosphere.  This  means  they  are  unstable  in  the  sense  used 
in  this  thesis.  That  physical  properties  of  the  lithosphere 
vary  sharply  has  been  inferred  from  the  geophysical 
observations.  If  the  plates  also  collide  with  each  other, 
inter-plate  forces  generate  concentrations  of  stress  on 
their  boundaries. 

Cox  (1973)  points  out  that  plate  tectonics  is  a 
unifying  hypothesis  which  provides  a  kinematic  model  of  the 
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upper  layer  of  the  earth,  and  that  it  can  be  used  to  make 
quantitative  predictions  about  most  phenomena  studied  by  the 
different  disciplines  of  the  earth  sciences.  Plate  tectonics 
links  all  these  phenomena  by  postulating  the  existence  of 
moving  plates.  Where  two  plates  are  moving  apart,  new 
oceanic  floor  is  formed  by  solidification  of  molten  rock  in 
the  opening  crack.  Where  two  converge,  one  of  them  is 
usually  thrust  beneath  the  other,  forming  oceanic  trenches, 
generating  earthquakes  and  as  the  descending  plate  melts 
rise  to  form  volcanoes. 

The  basis  of  this  hypothesis  has  been  developed  since 
the  beginning  of  the  century.  It  integrates  the  idea  of 
continental  drift  and  the  idea  of  sea-floor  spreading  (Hess, 
1962).  However,  the  hypothesis  itself  was  defined  during  the 
late  nineteen  sixties  (Le  Pichon,  1973)  when  the 
implications  of  the  following  ideas  were  understood: 

1.  The  primary  rheological  stratification  of  the  upper 
mantle  and  the  crust  into  lithosphere  and  asthenosphere 
governs  the  mechanical  behaviour  of  the  upper  layers  of 
the  earth  (Isacks  et.  al.,1968). 

2.  Most  of  the  mechanical  energy  of  the  surface  of  the 
earth  is  spent  within  a  few  narrow  seismic  belts 
(Richter ,  1 958 ) . 

3.  There  are  geometrical  constraints  imposed  on  the 
displacements  of  rigid  bodies  at  the  surface  of  the 
earth  (Wilson,  1965  and  Bullard,  1965). 
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In  the  development  of  the  plate  tectonics  hypothesis, 
there  are  many  other  important  contributions,  however,  it  is 
not  the  purpose  of  this  work  to  present  a  complete  review  of 
how  it  has  reached  its  present  form,  nor  to  go  into  the 
detailed  description  of  all  the  phenomena  that  contribute  to 
it.  Nevertheless,  the  aspects  mentioned  above  are  of  vital 
importance  to  what  will  be  treated  in  this  thesis.  In 
particular,  the  rheological  stratification  of  the  mantle  and 
lithosphere,  and  the  different  mechanical  properties  of  the 
diverse  intraplate  boundaries  are  important. 

The  following  concepts  must  be  underlined: 

1.  The  1 i thosphere , which  includes  the  crust  and  the 
uppermost  mantle,  has  significant  strength  and  may  be  up 
to  several  hundred  kilometers  in  thickness  (Isacks  et 
al.,  1968). 

2.  The  asthenosphere  extends  from  the  base  of  the 
lithosphere  to  a  depth  of  a  few  hundred  kilometers  and 
is  a  layer  of  effectively  little  strength  (Isacks  et 
al.,  1968). 

3.  A  transfom  fault  is  the  interplate  boundary  where  the 
relative  velocity  of  one  plate  with  respect  to  another 
is  parallel  to  the  boundary  itself  (Wilson,  1965). 
Turcotte  (1974),  explained  this  type  of  faults  as  the 
result  of  thermal  contraction  of  the  earth’s  crust,  and 
suggest  that  they  might  be  associated  with  the  formation 
of  a  solid  layer  due  to  the  cooling  of  a  nonelastic 
’’fluid".  Another  possible  explanation  is  that  transform 
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faults  are  a  product  of  intraplate  geometrical  evolution 
(McKenzie  and  Morgan,  1969). 

4.  Converging  zones  are  the  intraplate  boundaries  where  two 
plates  collide  with  each  other.  If  two  plates  of  the 
same  constitution  meet,  they  deform  and  can  produce 
mountain  belts,  like  the  Himalayas  produced  by  the 
collision  between  India  and  Asia. 

If  an  oceanic  plate  collides  with  a  continental  one 
then  the  heavier  oceanic  plate  subducts  beneath  the 
lighter  continental  one.  This  subduction  zone  is  a 
consuming  plate  margin  called  a  trench  (i.e.  McKenzie 
and  Morgan,  1969),  or  an  arc  (i.e.  Wilson,  1965)  which 
are  defined  as  lines  of  relative  motion  along  which 
surface  is  destroyed  asymmetrically.  The  surface  is 
destroyed  only  on  one  side  of  the  line. 

5.  Diverging  zones  or  mid-ocean  ridge  crests,  are  lines  of 
relative  motion  along  which  surface  is  produced 
symmetrically.  The  actual  relative  direction  of  motion 
need  not  be  perpendicular  to  this  line,  but  it  often  is 
(Le  Pichon ,1973). 

6.  Triple  junctions  are  points  on  the  surface  of  the  earth 
where  three  plates  meet.  The  relations  between  the 
relative  velocities  of  the  plates  determine  the 
evolution  of  the  triple  junction  in  time.  If  the 
evolution  of  the  triple  junction  is  possible  without  a 
change  in  geometry,  then  the  vertex  is  defined  as  a 
stable  junction.  If  the  junction  can  exist  for  an 
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instant  only  it  is  defined  as  unstable  (McKenzie  and 
Morgan,  1969).  The  junction  of  three  transform  faults  is 
always  unstable  whereas  the  junction  of  three  ridges 
spreading  perpendicularly  to  their  axis  is  always 
stable.  Junctions  with  two  boundaries  in  a  straight 
line,  fixed  with  respect  to  the  plate  they  bound,  are 
also  always  stable  (Le  Pichon,  1973).  An  unstable  triple 
junction  evolves  to  a  more  stable  kind  of  triple 
junction . 

It  is  also  necessary  to  differentiate  between  oceanic 
and  continental  lithosphere  because  their  material 
properties  are  different.  Most  investigators  who  have  tried 
to  define  the  detailed  mechanical  properties  of  the 
lithosphere  have  converged  to  a  model  of  oceanic  lithosphere 
that  is  consistent  with  the  thermal  models  of  the  upper 
mantle,  experimental  rock  mechanics,  gravity  and  bathymetry 
surveys,  and  patterns  of  seismicity  and  focal  mechanisms.  In 
this  model  the  lower  oceanic  lithosphere  is  visco-elastic 
and  will  not  support  loads  for  long  time  periods,  but  its 
stiffness  may  be  important  for  transient  loads  or  in  the 
initial  response  to  loading. 

The  upper  lithosphere  is  elastic  and  is  able  to  support 
loads  for  10’  years  without  appreciable  relaxation  (Forsyth, 
1979).  At  a  high  state  of  stress  it  deforms  elastically, 
yielding  brittle  fracture.  The  transition  from  the  brittle 
regime  of  the  upper  lithosphere  to  the  ductile  regime  of  the 
lower  layer  is  predominantly  thermally  controlled,  thus  the 
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flexural  rigidity  increases  as  the  lithosphere  cools.  This 
means  that  near  a  trench  brittle  failure  may  occur  at 
greater  depths  than  is  possible  near  the  ridges. 

This  pattern  of  higher  rigidity  for  the  cooler 

* 

lithosphere  would  imply  that  for  continental  lithosphere, 
which  is  older,  we  should  expect  higher  rigidity  values  than 
those  for  ocean  basins  provided  temperature  is  the  dominant 
parameter.  Nevertheless,  the  flexural  rigidity  values 
obtained  for  samples  of  continental  materials  are  smaller 
than  those  from  oceanic  lithosphere.  The  explanation  (e.g. 
For syth ,  1 979 ,*  McNutt,  1  980)  lies  in  the  mineralogic 
composition  of  the  continental  rocks  which  in  general  are 
weaker  and  subjected  to  continuous  isostatic  movements. 
Therefore,  the  ductile  regime  in  continental  areas  may  in 
general  be  reached  at  a  shallower  depth  than  in  the  oceanic 
crust . 

Observations  of  the  dispersion  of  Rayleigh  waves  have 
shown  a  significant  heterogeneity  in  the  structure  of  the 
upper  mantle  beneath  the  Pacific  Ocean.  The  heterogeneity  is 
well  correlated  with  lithospheric  age.  The  lid  to  the 
low-velocity  zone  is  very  thin  near  the  ridge  crest  and 
becomes  thicker  as  it  moves  away  from  the  ridge,  reaching  a 
thickness  of  about  60  km  after  50  million  years  and  about 
85  km  after  100  million  years  (Leeds  et.  al.,  1974). 

Hanks  and  Rayleigh  (1980)  introduce  the  notion  of 
differentiation  between  the  lithosphere  and  the  "elastic 
plate",  the  upper  part  of  the  lithosphere.  They  point  out 
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that  away  from  the  major  plate  boundaries,  the  elastic  plate 
can  support  several  kilobars  of  deviatoric  stress,  in 
response  to  local  or  regional  loads.  Thatcher  et  al.  (1980) 
by  looking  at  the  time-dependent  aseismic  deformation  of  the 
lithosphere  resulting  from  large  earthquakes  suggest  a 
thickness  of  the  elastic  plate  of  30  to  60  km;  which  is  done 
by  modelling  the  lithosphere  as  an  elastic  plate  overlying  a 
viscoelastic  half  space. 

McNutt  (1980)  concluded  from  gravity  observations  that 
the  upper  several  tens  of  kilometers  of  both  oceanic  and 
continental  intraplate  regions  can  support  a  kilobar  or  so 
of  deviatoric  stress  for  roughly  109  years  for  intact  rock. 
Turcotte  et  al.  (1980)  suggest  that  there  is  no  significant 
difference  in  elastic  strength  between  the  elastic  plates  of 
both  continental  and  oceanic  lithosphere,  at  least  in  the 
orogenically  stable  areas.  Moreover,  in  both  areas  the 
elastic  plate  has  significant  strength  and  a  thickness  of 
several  tens  of  kilometers,  on  a  time  scale  of  billions  of 
years . 

The  thickness  of  the  lithosphere  derived  from  flexural 
rigidity  studies  is  the  equivalent  thickness  of  a  plate  with 
uniform  properties,  responding  to  a  load  applied  on  it 
(Cazenave  et  al.,  1980).  The  values  derived  in  this  way  for 
the  thickness  of  the  oceanic  lithosphere  are  considerably 
less  than  those  determined  through  seismic  wave  propagation 
studies  (Leeds  et  al.,  1974).  A  comparison  between  elastic 
and  seismic  thickness  values  for  lithosphere  of  the  same  age 
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suggest  a  ratio  of  about  1/3  (Cazenave  et  al.f  1980). 

There  are  many  estimates  of  lithosphere  thickness  and 
they  need  not  be  the  same.  Seismic  lithosphere  is  defined  as 
the  layer  of  high  velocities  and  low  attenuation  of  seismic 
waves.  Another  definition  of  lithosphere  thickness  is  the 
thermal  thickness,  which  separates  the  region  where  the 
conductive  heat  transport  dominates  from  lower  regions  where 
convection  is  the  principal  mechanism  of  heat  transport 
(Parsons  and  Sclater,  1977). 

Based  on  the  existence  of  the  elastic  plate  described 
earlier,  the  determination  of  stresses  in  the  upper  few  tens 
of  kilometers  of  the  earth  can  be  approached  by  analysing 
elastic  models,  if  the  region  under  study  is  away  from  the 
interplate  boundaries  (Solomon  et  al.,  1980).  The  modelling 
of  regions  within  these  boundaries  should  be  done  by  taking 
into  account  the  thermal  dependence  of  the  stresses  that  may 
be  present  in  these  areas.  This  approach  will  be  followed 
for  the  different  geodynamic  problems  in  the  next  chapters. 
Kirby  (1977)  looks  into  this  problem  and  gives  the  following 
suggestions: 

1.  A  significant  portion  of  the  lithosphere  could  be 
described  using  a  nonlinear  strain-stress  relation.  At 
greater  depths  e  3  . 

2.  The  upper  part  of  the  lithosphere  should  deform  in  the 
brittle  failure.  He  mentions  that  in  particular,  in  a 
subduction  zone  the  upper  25  km  of  the  plate  should  be 
under  brittle  deformation. 
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3.  The  stresses  at  the  base  of  the  lithosphere  have  a  large 
range  of  variation  due  to  change  in  the  thickness  of  the 
plates . 

It  is  also  of  prime  importance  to  include  the  presence 
of  fluids  in  the  materials  modelled,  since  it  has  been  long 
recognized  that  variations  in  the  pressure  of  these  internal 
fluids  can,  and  in  fact  do,  change  the  state  of  stresses  in 
a  given  region.  This  has  been  observed  primarily  in  the 
induced  seismic  activity  that  has  occurred  in  some 
artificial  lakes  after  the  impounding  of  the  water  reservoir 
had  started  as  shown  by  Gough  (1978),  Simpson  (1976,  1981), 
Withers  and  Nyland  (1976,  1978),  Bell  and  Nur  (1978)  and 
many  others.  The  importance  of  this  pore  pressure  is 
dependent  on  the  porosity  and  permeability  of  the  rocks  in 
each  model. 

1.5  Methods  For  Determining  Tectonic  Stress  States. 

The  seismic  behaviour  of  different  geodynamic  systems 
varies  widely  from  relatively  small  microearthquakes  only 
detectable  with  highly  sensitive  instruments  to  large 
destructive  earthquakes.  The  kind  of  seismic  manifestations 
observed  in  a  particular  region  must  be  representative  of 
the  state  of  stress  that  prevails  in  it.  Different  seismic 
behaviour  should  then  characterize  locations  under  different 
stress  distributions. 

For  this  reason  I  have  started  this  thesis  with  a 
description  of  those  aspects  of  earthquake  seismology  that 
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will  allow  the  definition  of  a  stability  criterion 

representative  of  the  risk  of  an  earthquake.  This  criterion 

should  be  such  that  it  can  be  modified  to  become 

representative  of  the  different  stress  enviroments  which 

\ 

occur  in  the  upper  part  of  the  lithosphere. 

The  use  of  materials  with  a  linear  elastic 
stress-strain  relationship  will  prove  to  be  a  good  first 
approximation  in  the  selected  cases.  This  will  help  in  the 
treatment  of  more  complicated  situations.  Elastic  analysis 
combined  with  the  hydrological  characteristics  of  the  site 
can  result  in  reasonable  conclusions  about  the  distribution 
of  risk  in  some  areas.  The  modelling  of  materials  with 
thermoelastic  stress-strain  relationships  will  enable  us  to 
get  an  insight  into  the  physics  behind  the  observed  long 
range  precursors  of  some  large  earthquakes.  The  use  of  a 
porous  media  will  strongly  suggest  that  the  geometry  of  an 
applied  water  load,  its  loading  history  and  the  presence  of 
weakness  planes  are  the  parameters  that  very  well  determine 
the  presence  of  induced  seismicity  in  particular  areas.  This 
last  exercise  will  set  the  background  necessary  to  analyse 
the  classical  approach  to  consolidation  problems  in  water 
reservoir  engineering  in  order  to  seek  the  time  delays  which 
are  often  observed  in  induced  seismicity  cases.  The  usual 
analysis  does  not  predict  these  time  delays  in  a 
satisfactory  way. 

In  order  to  accomplish  any  of  these  goals,  it  is 
necessary  to  estimate  the  stresses  at  the  location  of 
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interest.  A  logical  way  to  proceed  is  to  try  to  make  some 
sense  of  the  state  of  stress  from  observed  measurable 
phenomena;  this  is  the  case  of  the  seismic  behaviour  which 
can  be  quantified. 

Two  admittedly  imperfect  methods  exist  for  assessing 
some  but  not  all  aspects  of  seismic  risk,  by  microseismicity 
studies  and  by  inference  of  the  stress  states.  The 
occurrence  of  microseismicity  reflects  stress  release  on  the 
faults  on  which  it  occurs.  The  frequency  and  size  of  such 
activity  must  depend  to  a  large  extend  on  available  energy. 
However,  the  ambient  stress  level  after  the  occurrence  of  an 
earthquake,  which  can  exist  without  being  reflected  in 
seismic  activity,  remains  unknown.  This  level  is  controlled 
by  fault  strength  and  can  vary  with  geology,  and  external 
influences  such  as  dams. 

The  magnitude-frequency  distribution  of  most  earthquake 
sequences  follows  the  empirical  rule 

Log  N=a-bM  (1.6) 

where  N  is  the  number  of  shocks  of  magnitude  greater  than  or 
equal  to  M.  "a"  and  "b"  are  constants  determined  for  each 
region.  The  constant  a  is  uncertain  since  it  depends  on  the 
sampling  time  and  is  subject  to  large  errors.  The  value 
usually  used  to  represent  the  state  of  stress  of  the  region 
is  "b".  Values  of  b  have  been  reported  from  0.5  to  1.5,  but 
usually  lie  between  0.7  and  1.0  for  tectonic  regions. 
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If  the  slope  of  equation  1.6  is  large  (b  large)  it 
predicts  a  smaller  ratio  of  large  to  small  earthquakes.  This 
suggests  the  presence  of  many  cracks  and  a  weak  material. 

Due  to  the  logarithmic  nature  of  equation  1.6  a  small  b 
value,  representative  of  a  zone  with  fewer  small  earthquakes 
but  perhaps  more  large  ones  will  indicate  a  stronger  (more 
homogeneous)  material,  since  in  general  there  are  fewer 
large  seismic  events  than  small  ones.  The  b  values  for 
aftershock  sequences  at  lakes  have  been  observed  to  be  about 
1.0,  and  they  are  much  higher  than  the  normal  b  values  of 
the  same  region.  The  b  values  of  the  foreshocks  are  even 
larger  for  large  induced  events. 

The  b  value  of  earthquakes  has  been  extensively  studied 

(Gupta  and  Rastogi , 1 976 ;  Gough, 1978)  and  has  been  found  to 

be  characteristic  of  the  seismicity  within  a  given  region.  A 

region  of  high  strength  and  variable  stress  distribution  is 

often  characterized  by  low  values  of  the  slope  b,  whereas  a 

source  region  in  which  there  are  existing  fractures  near 

* 

critical  stress  has  high  values  of  b  (Gough,  1978) 

Scholz  (1968)  relates  high  b-values  to  zones  under 
spatially  heterogeneous  stress  distributions  and  lower 
b-values  to  more  homogeneous  distributions.  He  concludes 
that  the  state  of  stress  plays  the  most  important  role  in 
determining  the  value  of  b.  This  means  that  for  an  area  with 
high  tectonic  activity  we  are  to  expect  a  high  b-value. 

Gupta  and  Rastogi  (1976)  and  Simpson  (1976)  explain  that  the 
high  b-values  observed  at  lake  sites  are  due  to  the 


.  ? ' 


25 


heterogeneous  nature  of  the  stresses  induced  by  the 
reservoir  depending  on  the  geological  structures  present. 

The  major  limitation  of  empirical  seismological 
relations  derived  from  the  study  of  seimic  catalogues  is 
that  these  are  not  by  any  means  complete;  that  is  the 
recording  network  is  only  capable  of  registering  all  the 
events  of  magnitude  greater  than  a  given  threshold.  Even  in 
the  best  studied  areas  the  period  when  records  are  available 
is  a  small  sample  in  geological  time  scales.  Nevertheless, 
that  is  all  we  have  to  go  on.  Lamoreaux  et  al.  (1983) 
studied  this  problem  of  incomplete  data  in  their  study  of 
the  occurrence  of  seismic  patterns  in  the  Middle-Amer ica 
trench.  Discussions  about  the  validity  of  statistical 
seismic  studies  have  been  done  by  many  independent 
researchers  (ie.  Lomnitz,  1966;  Lamoreaux,  1982;  Rikitake, 
1977).  They  all  concluded  that  due  to  the  short  history  of 
recording,  statistics  may  be  misleading,  nevertheless,  that 
is  all  we  have  to  go  on.  Rikitake  (1977),  points  out  that  it 
is  possible  to  use  statistics  in  seismology  only  if  a  series 
of  seismic  observations  over  a  period  of  several  decades  or 
longer  is  available  for  that  particular  region. 

Geologic  investigation  in  the  late  Quaternary  might 
provide  valuable  data  on  recurrence  periods  of  damaging 
earthquakes  (Allen,  1975).  Such  investigations  appear  to 
help  overcome  many  of  the  statistical  shortcomings  of 
instrumental  and  historical  records.  It  may  be  expected  that 
the  dangerous  large  events  (M>7)  will  leave  geologic  traces, 
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from  which  not  only  their  occurrence,  but  also,  the  amount 
of  slip  and  thus  their  magnitude  can  be  derived.  These 
investigations  may  provide  the  "maximum  credible  earthquake" 
data  required  for  hydroelectric  project  design,  for  other 
large  engineering  projects  and  is  desired  in  the  assesment 
of  earthquake  risk  in  geodynamical  modelling. 

A  third  method  of  stress  analysis  is  to  attempt  a 
direct  measurement  of  stress.  Several  methods  exist  for 
routine  determination  of  formation  stress  (Haimson  and 
Fairhurst,  1970).  The  most  useful  appears  to  be  that  used  by 
the  U.S.  Bureau  of  Mines  (Obert  et  al.,  1962).  This  method 
involves  placing  a  diameter  gauge  in  a  borehole  of 
approximately  3  cm  diameter.  This  is  usually  a  low  modulus 
inclusion  which  consists  of  a  piston  in  a  brass  barrel.  The 
piston  measures  the  borehole  diameter  by  deflection  of  a  bar 
fitted  with  a  strain-gauge.  The  resistance  of  the  gauge  can 
be  monitored  with  normal  strain  gauge  equipment. 

Once  the  borehole  gauge  is  in  place  a  large  diameter 
coring  bit  20  cm)  can  be  used  to  drill  out  the  rock 
containing  the  gauge.  A  fairly  simple  application  of 
elasticity  theory  relates  the  observed  relaxation  of  the 
rock  to  stresses  in  the  plane  perpendicular  to  the  borehole. 
Three  measurements  of  this  type  in  3  independent  directions 
are  used  to  deduce  the  complete  stress  tensor.  All  that  is 
required  is  a  competent  rock  drilling  crew  with  equipment 
capable  of  cutting  large  diameter  cores.  It  is  crucial  that 
the  drilling  crew  be  experienced  for  the  large  diameter  core 
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must  be  retrieved  in  one  piece.  Borehole  gauges  can  be 
constructed  by  any  qualified  machinist  and  strain  monitoring 
equipment  is  a  standard  item.  The  problem  is  that  these 
measurements  are  local  and  shallow,  and  not  necessarily 
representative  on  a  regional  scale.  The  seismic  behaviour  on 
the  other  hand  is  regional  in  character  and  can  be  thought 
as  a  way  of  determining  the  average  behaviour  of  a  given 
zone . 


Flat  jack  tests  are  not  as  reliable  for  these  stress 
analyses.  The  results  can  be  severely  contaminated  by 
geometric  influences  due  to  the  excavation  from  which  they 
are  made.  High  modulus  inclusion  methods  (Hast,  1958)  might 
have  a  place  in  competent  rock  but  seem  ruled  out  in 
sedimentary  basins  or  heavily  fractured  materials. 

Hydrofracturing  (Hubbert  and  Willis,  1957;  Haimson  and 
Fairhurst,  1970)  is  probably  the  only  way  to  acquire  stress 
data  at  depth  but  it  requires  the  drilling  of  several  deep 
wells.  Such  wells  would  however  serve  a  useful  purpose  in 
monitoring  fluid  pressure  behaviour  during  loading  of  a 
reservoir.  This  technique  has  limitations 

1.  One  principal  stress  must  be  vertical. 

2.  Drill  hole  must  be  near  vertical. 

3.  Cracks  can  introduce  considerable  uncertainty  in 
results . 

The  orientation  of  the  principal  stresses  in  the  upper 
few  kilometers  of  the  earth's  crust  has  also  been  determined 
by  analysing  the  direction  of  oil-well  breakouts  (i.e.  Gough 
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and  Bell,  1981).  This  is  done  with  a  4-arm  dipmeter,  which 
records  the  magnitude  and  azimuth  of  two  orthogonal 
diameters  of  the  hole.  A  breakout  occurs  when  when  one  of 
the  diameters  fits  the  original  drill  and  the  other  is 
greater.  The  observed  direction  of  the  breakouts  is 
consistent  within  a  single  well,  and  seems  to  be  independent 
of  the  geological  structures  through  which  the  hole  has  been 
drilled.  The  consistency  of  breakout  orientation  throughout 
a  single  geologic  province  (Bell  and  Gough,  1979),  suggests 
that  the  systematically  oriented  azimuths  are  products  of 
stress  concentrations  near  the  hole  in  a  regional  stress 
field  where  the  two  horizontally  oriented  principal  stresses 
are  different . 

Well  breakouts  give  no  measure  of  the  magnitude  of  the 
stresses  and  if  they  are  to  be  useful  in  the  determination 
of  the  regional  orientation  of  the  principal  stresses  we 
require  the  presence  of  many  wells  in  the  area. 

Finally,  one  of  the  methods  for  estimating  stresses  in 
the  lithosphere  is  by  modelling,  using  the  major  geophysical 
observables  of  heat  flow,  seismic  wave  propagation  data, 
gravity,  radioactivity,  electric  conductivity,  magnetometry , 
and  the  plate  velocities  as  constraints.  The  link  between 
the  geophysical  observables  and  the  stress  in  the 
lithosphere  is  through  the  different  constitutive  equations 
mentioned  earlier  in  this  chapter.  These  are  equations  that 
relate  the  rate  of  strain  to  the  applied  stress, 
temperature,  pressure  or  any  other  parameter  that  might 
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influence  this  strain  rate.  Unfortunately  almost  all  these 
parameters  must  be  estimated.  Model  calculations  are 
discussed  later  in  this  thesis. 

1.6  Some  Earthquake  Related  Phenomena 

Several  aspects  of  the  earth's  seismic  behaviour,  that 
are  of  importance  to  the  understanding  of  this  work  and  are 
the  necessary  background  to  the  following  chapters  must  be 
outlined . 

Recurrence  patterns. 

The  study  of  the  seismic  catalogs  for  the  different 
parts  of  the  earth  has  resulted  in  several  empirically 
derived  observations.  One  of  these  is  that  large  destructive 
earthquakes  tend  to  repeat  themselves  in  location  and  in 
time  (i.e.  McNally  and  Minster,  1981;  Singh  et  al.,  1981). 
For  simple  plate  boundaries  the  average  recurrence  time  7 
can  be  calculated  from 

T=u/av  (1.7) 

where  u  is  the  average  displacement  that  takes  place  at  a 
location  due  to  seismic  activity,  a  is  the  ratio  of  seismic 
slip  to  the  total  slip  and  v  is  the  relative  plate  velocity 
(Sykes  and  Quittmeyer,  1981). 

Two  assumptions  are  implicit  in  this  kind  of  formulation, 
one  is  that  the  major  part  of  the  slip  at  the  location  is 
due  to  large  earthquakes  and  the  other  is  that  the  plate 
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velocities  can  be  treated  as  constant  over  large  periods  of 
time  (Sykes  and  Qui ttmeyer , 1 98 1 ) .  This  has  lead  to  the 
development  of  three  different  type  of  recurrence  patterns 
(Shimazaki  and  Nakata,  1980).  These  patterns  (figure  2)  are: 

1.  The  time  recurrence  model  shown  in  figure  2. a  assumes 
that  the  stress  required  to  fracture  a  given  location  is 
constant  as  indicated  by  the  horizontal  line  at  Tt  level 
of  stress  on  this  figure.  This  is  not  necessarily  true 
since  fractures  decrease  the  strength  of  the  material 
and  modifies  its  frictional  resistance.  The  model 
implies  that  stress  is  accumulated  until  it  reaches 
level  T,.  Then  an  earthquake  takes  place  and  the  stress 
level  drops  to  a  smaller  level  and  the  cycle  starts 
again.  This  pattern  enables  us,  in  principle,  to  predict 
when  is  an  earthquake  going  to  occur  based  on  the  size 
of  the  previous  one. 

2.  Figure  2.b  describes  the  Slip  predictable  model.  It 
assumes  that  the  level  of  stress  after  the  occurrence  of 
any  seismic  event  is  the  same,  that  if  an  earthquake  is 
going  to  occur  today,  its  magnitude  can  be  predicted 
from  the  time  lag  from  the  previous  one.  This  kind  of 
modelling  is  not  in  good  agreement  with  the  seismic 
catalogues  from  Japan  (Shimazaki  and  Nakata,  1980) 
California  (Bufe,  1977)  and  the  Middle  America  Trench 
(Singh  et  al . ,  1983). 

3.  The  strictly  repetitive  model  shown  in  figure  2.c  is  the 
combination  of  the  previous  two,  the  levels  T,  and  T2 


'  A. 


STRESS 


31 


a 
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Earthquake  recurrence  patterns 
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are  fixed,  the  time  of  recurrence  is  always  the  same  and 
the  earthquake  is  of  the  same  magnitude  every  time.  This 
one  has  all  the  defects  of  the  other  two  combined,  and 
is  considered  a  limiting  ideal  case. 

The  common  assumption  of  all  these  models,  that  the 
rate  of  accumulation  of  stress  is  constant,  has  been  used 
extensively  in  many  other  geophysical  studies  (i.e.  Newman 
and  Knopoff,  1982;  Singh  et  al.,1983).  This  concept  of  a 
continually  increasing  stress  will  permit  the  study  of  the 
variation  of  the  stability  function  due  to  an  earthquake 
(Chapter  3).  Therefore,  I  have  assumed  that  earthquakes 
represent  a  sharp  reduction  in  stress  and  are  the  only 
natural  way  in  which  stress  can  change  abruptly. 

Observed  recurrence  times  at  a  given  location  are 
irregular  enough  that  for  most  regions  the  knowledge  of  this 
average  time  interval  between  succesive  earthquakes  does  not 
greatly  constrain  future  times  of  occurrence  (Rikitake, 
1977). 

Se  i s/77  / c  / ty  Patterns . 

Another  observation  from  historical  seismicity  records 
is  that  sometimes  large  earthquakes  appear  to  affect  the 
seismic  behaviour  of  zones  which  are  very  far  away,  and  some 
large  events  seem  to  be  preceded  by  anomalous  seismic 
behavior.  All  these  types  of  peculiar  seismic  behaviour  have 
been  grouped  as  Seismicity  Patterns  or  Long  Range 
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Most  of  the  observed  seismicity  patterns  can  be 
interpreted  as  a  combination  of  quiescence  and  activation 
(Lamoreaux,  1982).  Quiescence  is  understood  as  a  decrease  in 
the  seismic  activity  of  a  region  in  comparison  to  normal 
activity.  Activation  is  then  an  increase  in  the  seismicity 
of  a  region,  an  increase  in  the  total  energy  release  rate 
there,  or  both. 

Many  of  these  patterns  have  been  identified,  some  of 
them  involve  activation  only.  Foreshocks  are  the  most 
obvious  of  this  kind  of  precursors,  they  are  swarms  or  small 
events  close  in  time  and  space  that  occur  just  before  a 
major  earthquake  takes  place  in  the  region.  Successful 
predictions  have  been  made  by  recognizing  this  pattern  i.e. 
the  Haicheng  magnitude  7.3  earthquake  of  1975  (Scholz, 

1977).  Pattern  'swarms'  (Keilis-Borok  and  Rotvain,  1979) 
refers  to  a  group  of  small  magnitude  events  closly  grouped 
in  space  and  time  occurring  during  a  time  interval  when  the 
overall  seismicity  is  not  below  average,  and  takes  place 
several  years  before  a  major  earthquake  takes  place  in  the 
area.  McNally  (1977)  observed  this  pattern  in  California, 
Ohtake  (1976)  reported  similar  findings  in  Japan. 

Keilis-Borok  et  al.  (1980  a,b,c)  have  proposed  three 
related  seismicity  patterns  which  deal  with  the  activation 
of  a  region.  These  three  patterns  have  been  collectively 
termed  "Bursts  of  seismicity".  They  consist  of  abnormal 
clustering  in  time,  energy  and  space  before  a  major 
earthquake.  These  are,  pattern  'swarms',  pattern  'bursts  of 
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aftershocks',  which  consist  of  an  anomalous  number  of 
aftershocks  concentrated  at  the  beginning  of  the  aftershock 
sequence,  and  pattern  'Sigma'  which  consists  of  an  increase 
in  the  cumulative  seismic  energy  released  to  the  2/3  power 
over  a  sliding  time  window,  pattern  sigma  is  identified  as 
the  peak  of  the  summation. 

Some  of  the  observed  seismicity  patterns  are 
superposition  of  quiescence  and  activation.  Evison  (1977), 
suggested  a  pattern  of  premonitory  swarms  followed  by  a 
quiescent  period  in  the  epicentral  region.  Mogi  (1969) 
identified  a  doughnut  pattern  for  some  large  earthquakes  in 
Japan.  Keilis-Borok  and  Rotvain  (1979),  suggest  a  similar 
algorithm,  pattern  'A-Q'  (activation-quiescence).  Pattern 
'A-Q'  combines  two  activation  areas  at  opposite  sides  of  an 
earthquake  prone  site  with  a  quiescent  area  in  the  middle. 

There  is  no  satisfactory  physical  explanation  for  the 
occurrence  of  these  patterns,  although  some  theories  exist. 
Two  major  approaches  have  been  developed 

1.  They  are  a  result  of  the  actual  distribution  of  physical 
properties  of  the  earth. 

2.  They  are  the  result  of  the  processes  of  crack  population 
growth,  independent  of  what  material  these  cracks  are 

on . 

The  right  answer  will  probably  link  both  of  these 
explanations.  In  chapter  4  this  is  explored  by  using  the 
stability  functions. 


2.  Types  of  failure  and  stability  studies. 

A  major  aim  of  this  work  is  to  link  together  what  is 
known  about  the  earth  processes  from  direct  in-situ 
measurements  and  experimental  results  with  the  variation  of 
the  geophysical  observables  such  as  se i smological  behaviour. 

To  accomplish  this  it  is  convenient  to  review  some 
concepts  of  rock  mechanics  which  are  fundamental  to  the 
further  understanding  of  what  is  expected  from  the  earth. 
This  analysis  will  set  up  the  rules  that  must  be  followed  in 
any  kind  of  description  of  rock  failure. 

Once  the  mechanical  background  has  been  presented,  the 
chapter  proceeds  to  explore  different  ways  to  determine 
seismic  risk.  Applications  to  time  dependent  phenomena  are 
analysed  first.  This  is  modified  for  the  determination  of 
seismic  risk  isolines  in  active  areas.  The  chapter  concludes 
with  a  method  in  which  the  seismic  risk  is  represented  by 
probability  of  seismic  failure  in  particular  areas. 

2 . 1  Rock  Failure . 

It  is  important  to  recall  the  different  types  of 
failure  criteria  that  are  applied  in  rock  mechanics.  All 
these  are  functions  of  the  internal  distribution  of  stress 
in  the  material.  However,  the  internal  stress  of  small 
laboratory  samples  may  not  be  representative  of  the  stress 
distribution  in  real  geological  environments.  Nevertheless, 
all  the  criteria  assume  that  the  "average"  stress  is 
representative  and  can  be  used  to  describe  this  behaviour. 
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Von  Mises  criterion  suggests  that  yielding  of  ductile 
material  occurs  when  the  limit  value  of  the  octahedral  shear 
stress  ( r © )  is  achieved  (Jaeger  and  Cook,  1979).  r0  as  a 
function  of  the  principal  stresses  alf  o2  and  o3  is  given 
as: 

T0=‘\/3{(o)-o2)2  +  (o2-a3)2  +  (o3-o1)2}l/2  (2.1) 

Another  yielding  criterion  is  that  due  to  Tresca,  which 
states  that  ductile  yielding  will  occur  when  the  maximum 
shear  stress  rmx  is  reached  (Jaeger  and  Cook,  1  979),  where*. 

Tmx=(o ,-a 3) /2  (2.2) 


The  most  widely  used  but  not  necessarily  the  most 
appropriate  failure  criterion  is  the  one  known  as 
Mohr-Coulomb .  This  criterion  implies  that  failure  occurs 
when  the  critical  shear  stress  r,  given  as  r=So+tftan0,  is 
reached.  S0  and  0  are  material  constants.  This  criterion 
(figure  3),  will  be  discussed  extensively  in  the  following 
parts  of  this  chapter.  In  this  figure  are  shown  the 
classical  ways  in  which  this  criterion  can  be  represented, 
a)  In  shear  stress  vs  normal  stress  in  cartesian 
coordinates,  and  b)  in  a  principal  stress  diagram,  where 
m= ( 1 +sin0 ) /( 1 -sin0 )  and  b=2Socos0/( 1 -sin0 ) . 

In  three  dimensions  the  stress  envelope  of  figure  3.b 
is  given  by  the  hexagonal  pyramid  shown  in  3.c. 
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Figure  3....  Mohr-Coulomb  failure  criterion  in  a  shear 
stress  vs  normal  stress  diagram  (a)  and  in  principal 
stresses  diagrams  in  two  dimensions  (b)  and  in  three 
dimensions  (c ) . 
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The  Drucker-Prager  criterion  is  based  on  the 
replacement  of  the  pyramid  of  figure  3.c  from  the 
Mohr-Coulomb  method  by  a  circular  cone  inside  the  pyramid. 
This  results  in  the  following  failure  criterion 

F  =  aJ1+i/J2-k  =  0  (2.3) 

where  Ji  and  J2  are  the  first  and  second  stress  invariants 

a=tan0/(9+12tan20) 1 '  2 
and 

k=3So/(9+12tan20) 1 '  2 

The  importance  of  this  kind  of  criterion  (figure  4),  is 
that  it  takes  into  account  the  intermediate  principal  stress 
a2,  which  is  neglected  in  the  Mohr-Coulomb  criteria. 

All  these  criteria  have  been  postulated  on  the  basis  of 
theoretical  considerations  and  have  been  extensively  tested. 
Except  for  the  Mohr-Coulomb  their  application  in  geophysics 
is  rare  since  they  all  imply  assumptions  about  the  behaviour 
of  the  materials  of  interest.  A  failure  criterion  based  on 
the  actual  physical  process  leading  to  failure  was  suggested 
by  Griffith  in  1924.  A  material  contains  a  number  of 
randomly  oriented  zones  of  potential  failure  in  the  form  of 
weak  grain  boundaries,  fissures,  cracks  or  open  flaws.  These 
flaws  can  be  approximated  by  long  thin  elliptical  openings 
and  one  can  assume  that  there  is  no  interference  between 
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Figure  4 
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Drucker-Prager  failure  criterion. 
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flaws  in  an  idealized  elastic  material.  In  this  case  even 
under  high  compressive  states  of  stress,  tension  stresses 
occur  in  the  boundaries  of  these  openings  and  fracture  must 
initiate  at  that  boundary  when  the  tensile  strength  is 
exceeded . 

The  maximum  tangential  stress  (atn)  near  the  tip  of  the 
elliptical  opening  is 

o  t  n  *m=(7y±  (oy  2+t  x  y  2  )  1  '  2  (2.4) 

where  m=a/b. 

This  criterion  implies  several  characteristics  of  the 
crack  (Stagg  and  Zienkiewicz,  1968): 

1.  It  propagates  out  of  the  plane  of  the  flaw. 

2.  It  propagates  to  infinity  if  it  is  not  confined  (tensile 
splitting  rupture). 

3.  It  propagates  towards  the  major  principal  axis. 

4.  No  rupture  occurs  when  a,  and  a3  are  compressive  (stable 
crack  propagation). 

5.  The  final  length  of  the  crack  depends  on  the  ratio  o2/o^ 
and  the  initial  flaw  length. 

6.  Depending  on  the  crack  density,  cracks  may  merge  and  a 
failure  may  occur  due  to  this  crack  coalescence 
(creation  of  a  shear  fault). 

Hoek  and  Brown  (1981)  have  made  an  empirical  attempt  to 
put  all  the  information  available  together  (figure  5),  and 
develop  an  empirical  strength  criterion: 
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Figure  5....  Shows  some  examples  of  Hoek  and  Brown's 
empirical  failure  criterion  for  a)  sandstones  and  b)for 
granites  .(modified  from  Hoek  and  Brown  »  1981) 
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o i /o  0=o  3/o  o  + (Ma 3/a0  +  S ) 1 7 2  ^  (2.5) 

where  M  and  S  are  constants  which  depend  on  the  rock.  For 
intact  rock  S=1  and  M=a0/atf  o t  being  the  tensile  strength 
of  the  material  and  o0  the  confining  pressure.  M  will 
decrease  rapidly  as  the  degree  of  prior  fracturing 
increases . 

2.2  The  Stability  of  Earth  Materials. 

The  concept  of  stability  has  been  used  in  rock 
engineering  to  predict  failure.  In  major  engineering 
projects,  it  is  of  great  importance  to  determine  the 
strength  of  the  rocks  at  a  site  and  avoid  the  collapse  of 
the  geological  structures  that  may  be  present  there.  This  is 
done  by  making  several  assumptions  about  the  strength  of  the 
materials,  discontinuities,  and  other  local  features.  In 
general,  these  calculations  demonstrate  the  amount  of 
strengthening  required  in  the  area.  Introducing  artificial 
supports  will  locally  increase  the  material  strength  and 
permit  the  work  to  be  done  with  a  safety  factor. 

The  simplest  example,  involving  this  type  of  reasoning 
is  the  determination  of  the  average  cohesion  necessary  to 
maintain  the  equilibrium  of  an  open  pit  (Brown  and  King, 
1966).  They  determine  the  components  of  stress  and  make  use 
of  the  Coloumb  failure  envelope  to  determine  the  STABILITY 
of  a  slope  model,  by  looking  at  the  toe  of  the  slope,  and 
calculating  the  plane  of  possible  failure  by  using 
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8-  tan"12axy/  2(axx-ayy)  +tt/4  +0/2  (2.6) 

for  axx  greater  than  ayY  as  a  function  of  the  angle  of 
internal  friction  of  the  material  0. 

Once  the  potential  slip  surface  has  been  given  then  the 
average  cohesion  S0  for  each  path  L  can  be  calculated  ass 

S0=[ JrdL+J (a-P) tan0-dL]//dL  (2.7) 

where  r  is  the  shear  stress,  o  the  normal  stress  and  P  the 
pore  water  pressure.  The  Couloumb  failure  envelope  is  given 
as : 

| t  | >S0~ (a-P) tan0  (2.8) 

The  calculation  assumes  that  rock  behaviour  is 
linear-elastic  and  that  failure  will  occur  in  the  direction 
of  an  angle  8  for  which  the  corresponding  value  of  S0  is 
minimum.  This  engineering  approach  works.  In  this  thesis  we 
extend  the  analysis  to  the  study  of  stability  of  large 
geologic  structures  in  order  to  connect  seismicity  and  rock 
mechanics  in  a  quantitative  way. 

In  geodynamical  studies  of  the  elastic  plate, 
instabilities  represent  seismic  events  and  nothing  can  be 
done  to  increase  the  mechanical  strength  of  the  material. 
Turcotte  et  al.  (1980),  studied  instabilities  due  to  fault 
weakening  caused  by  frictional  work.  Stuart  (1979)  points 
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out  that  the  subject  of  earthquake  mechanics  viewed  as  a 
boundary  problem  in  continuum  mechanics  reduces  to  the  study 
of 

1.  constitutive  properties  and  geological  structures 
present  in  the  surroundings,  and 
2c  remotely  applied  boundary  conditions. 

These  concepts  must  lead  to  the  observed  nature  of  the 
earthquake  cycle.  The  stages  of  this  cycle  are 

1.  slow  increase  of  stress  due  to  remote  loading, 

2.  the  onset,  propagation  and  cessation  of  the  rupture  and 

3.  post  seismic  adjustment. 

To  this,  it  should  be  added,  that  fast  stress  changes  may 
occur  in  the  upper  few  kilometers  of  the  earth.  This 
inclusion  will  in  fact  integrate  the  engineering  approach 
and  the  more  general  geophysical  studies  of  the  earth's 
crust  and  upper  mantle. 

Earth  instability  models  require  a  precise  statement  of 
the  dependence  of  the  stability  on  the  material  properties 
of  the  models.  The  instability  criterion  depends  on  the 
stiffness  of  the  material  and  the  stiffness  of 
discontinuities  present  in  the  model.  It  is  given  as  the 
possibility  of  a  fault  slip  due  to  a  remote  displacement  or 
stress  release.  Stuart  (1979)  points  out  that  physically, 
this  condition  corresponds  to  a  sudden,  finite  fault  slip 
caused  by  an  arbitrarily  small  progressive  change  on 
regional  stress.  The  unstable  slip  during  faulting  is  driven 
by  the  release  of  elastic  energy  stored  in  the  surroundings 
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as  proposed  in  the  elastic  rebound  theory.  Then,  any 
constitutive  law  in  which  stress  can  decrease  rapidly  during 
deformation  can  be  used  to  create  an  instability  criterion. 


Failure  represents  the  actual  stress  release  and 
originates  at  the  location  in  which  the  ratio  of  the 
strength  of  the  material  to  the  stress  stored  in  it  is 
minimum.  The  crack  will  propagate  through  the  rock  until  it 
reaches  a  location  where  the  stresses  are  not  enough  to 
overcome  the  rock  strength.  Therefore,  the  size  of  the 
stress  concentration  will  determine  the  size  of  the  crack. 
Very  localized  concentrations  then  are  expected  to  produce 
small  cracks,  while  regional  concentrations  are  more  likely 
to  produce  large  cracks. 

The  stress  concentration  due  to  the  propagation  of  the 
seismic  waves  ahead  of  the  tip  of  a  crack  can  cause  the 
fault  slip.  The  cohesive  force  of  the  material  where  the 
rupture  occurs,  determines  what  happens  when  the  rupture 
propagates  along  a  fault  plane  with  obstacles  or  barriers. 
These  barriers  may  be  thought  of  as  localized  high  values  of 
the  cohesive  force  of  the  material.  Aki  and  Richards  (1980) 
suggest  that  three  different  things  may  occur  when  the  tip 
of  a  crack  passes  through  such  a  barrier,  depending  on  how 
strong  the  barrier  is  with  respect  to  the  initial  stress; 

1.  If  the  initial  stress  is  relatively  high,  the  barrier  is 
broken  immediately. 

2.  If  the  barrier  strength  is  relatively  high,  the  barrier 
is  left  behind  unbroken. 
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3.  If  they  both  are  similar,  the  barrier  is  not  broken  when 
the  crack  tip  passes,  but  it  may  eventually  break  due  to 
later  stress  increases. 


2.3  General  definition  of  stability. 

The  determination  of  zones  of  seismic  risk  in 
geosciences  is  a  subject  in  which  several  items  must  be 
taken  into  account: 

1 .  Concentration  of  stresses  due  to  geological  structures 
in  the  region  of  interest,  that  is,  sudden  changes  in 
material  properties  occurring  inside  the  model  under 
study . 

2.  External  forces  acting  on  the  medium,  like  plate  driving 
forces,  volcanic  activity,  igneous'  intrusions ,  changes 
in  water  content  or  in  pore  pressure  on  the  rock  mass. 

3.  Evidence  of  previous  seismic  activity  in  the  area. 

4.  Location  of  the  site  with  respect  to  the  active  seismic 
belts  of  the  earth. 

Once  all  the  above  factors  have  been  considered, 
assumptions  about  the  constitutive  equations  that  regulate 
the  stress  propagation  in  the  model  can  be  made.  The  next 
step  is  to  calculate  the  stress  tensor  components.  This  can 
be  done  by  solving  the  boundary  value  problem  with 
analytical  methods  or  numerical  techniques.  The  detailed 
development  examples  of  both  approaches  can  be  found  in 


appendices  A  and  B. 
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If  the  variation  of  the  stress  tensor  throughout  the 
model  is  known,  its  components  can  be  combined  in  such  a  way 
that  a  measurement  of  stability  at  each  point  of  the  model 
can  be  calculated.  This  is  done  by  applying  the  Mohr  circle 
principle  to  the  resulting  principal  stresses  and  relating 
it  to  a  failure  criterion. 

The  physical  significance  of  using  a  Coulomb  failure 
envelope  is  that  the  rock  is  taken  to  have  a  linear 
relationship  between  the  critical  shear  stress  r13  and  the 
critical  maximum  stress  diference  o ^-o 2  ,  that  is  required  to 
reach  failure.  A  solid  exhibiting  this  kind  of  behaviour  is 
characterized  by  a  linear  relation  between  the  components  of 
stress  and  strain  as  described  in  section  1.1  and  implies 
reversibility  of  the  stress  strain  relations.  This 
assumption  for  the  upper  few  kilometers  of  the  earth  is 
plausible  if  the  constraints  mentioned  through  chapter  1  are 
observed . 

If  the  envelope  of  the  Mohr  circles  at  failure  has  the 
form  of  a  straight  line  (figure  6),  it  is  given  by  equation 
(2.8),  which  can  be  rewritten  as: 

Ar=So+Aantan0  (2.9) 

where  A r,  is  the  shear  strength  increment  of  the  rock.  0,  is 
the  angle  of  shear  resistance,  A an,  is  the  normal  effective 
stress  increment  on  the  plane  of  fracture,  and  S0  is  the 
shear  strength  increment  of  the  material  under  zero  normal 
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Figure  6....  Shows  the  general  definition  of  instability  as 
the  minimum  distance  from  the  surface  of  the  Mohr  circle  to 
the  failure  criterion. 
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pressures . 

'S0’  varies  considerably  from  zero  in  fractured 
material  to  -  50  bars  for  sedimentary  rocks  and  up  to 
several  hundreds  of  bars  for  igneous  rocks  in  intact 
material  (Withers,  1977).  In  the  case  of  our  geological 
situation,  we  have  the  presence  of  fractures  and  *S0'  is 
expected  to  be  small.  The  shear  angle  <p  lies  between  25°  and 
45°,  therefore  the  coefficient  of  friction  is  between  0.47 
and  1.0.  This  one  is  usually  around  0.6  (0=30°). 

Since  the  value  of  'S0?  is  unknown,  fixing  it  to  zero, 
and  determining  the  variation  of  the  minimum  distance 
between  this  failure  envelope  and  the  Mohr  circle  will 
result  in  a  measurement  of  the  stabilty  of  the  system  at  a 
given  time.  This  is  equivalent  to  the  introduction  of  a 
reference  frame.  That  is,  at  the  beginning  of  the 
calculation,  the  model  was  taken  in  equilibrium,  and  then  by 
making  some  physical  change,  increments  on  the  principal 
stresses  are  calculated.  These  will  result  in  a  unique 
measurement  that  will  denote  if  the  corresponding  location 
in  the  model  has  increased  or  decreased  in  stability  with 
respect  to  the  original  situation. 

This  kind  of  reference  frame  can  only  be  used  if  there 
is  a  physical  change  in  the  conditions  of  the  area  from  one 
time  to  another.  The  risk  of  a  seismic  event  should  depend 
on  the  minimum  distance  between  the  failure  envelope  and  the 
resulting  Mohr  circle.  Some  of  its  applications  will  be  seen 
in  the  succeding  chapters  of  this  thesis. 
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Although  the  definition  of  stability  functions  requires 
the  precise  definition  of  a  failure  criterion,  it  is  not 
restricted  to  the  use  of  a  linear  envelope.  In  order  to 
apply  this  kind  of  stability  analysis  to  greater  depths  than 
those  for  which  elastic  behaviour  is  justified  the  failure 
criterion  must  be  changed.  In  Chapter  4  a  different  type  of 
envelope  has  been  used  since  the  stress-strain  relations 
dealt  with  there  are  not  always  linear,  and  the  material  is 
allowed  to  become  somewhat  ductile.  Therefore  the  envelope 
there  will  describe  the  yielding  of  rocks. 

The  shapes  of  the  reference  failure  envelopes  used 
throughout  the  following  chapters  are  derived  from  the 
physical  characteristics  that  will  be  described  for  each 
model.  In  general  the  materials  are  assumed  to  have  zero 
tensile  strength,  to  behave  elastically  in  the  upper  few 
kilometers  of  the  earth's  lithosphere  and  the  shear  stress 
angle  <t>  is  taken  as  30°.  As  the  temperature  of  the  model 
reaches  300°C  plastic  behaviour  is  assumed. 

* 

2.4  Instability  function  in  three  dimensions. 

In  the  Mohr's  circle  representation  in  three 
dimensions,  the  normal  and  shear  stress  across  a  plane  of 
weakness  whose  normal  has  direction  cosines  1,  m,  n,  are 
given  by  Jaeger  and  Cook  (1979,  p27).  Fixing  two  of  the 
direction  cosines  (say  n  and  1)  two  equations  can  be 
obtained.  Each  of  them  represents  one  family  of  Mohr’s 
circles  in  two  dimensions  and  for  a  fixed  value  of  the 
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corresponding  direction  cosine  each  represents  a  unique 
circle.  Therefore,  by  fixing  n  and  1,  two  circles  can  be 
drawn  such  that  their  intersection  will  lie  at  a  point  on 
the  surface  of  a  three  dimensional  Mohr  representation,  and 
will  be  a  unique  location  for  these  two  circles  whose 
centers  are  (a1+a2)/2  and  ( o2+oz)/2  and  whose  radii  are  AC 
and  BD  respectively  (figure  7). 

With  the  previous  procedure  it  is  possible  to  determine 
the  values  of  o  and  r  for  every  combination  of  stresses. 

That  is,  the  location  of  the  point  P  from  figure  9  can  be 
determined  for  any  time  using  the  simple  Mohr-Coulomb 
failure  envelope.  This  criterion  suggests  that  failure 
occurs  when  the  minimum  shear  on  a  failure  plane  exceeds  the 
shear  strength  of  the  rock. 

If  fractures  are  present  in  the  material  S0  is  probably 
small.  Therefore  since  the  value  of  S0  is  unknown,  I  set  it 
to  zero.  Now  the  variation  of  the  minimum  distance  between 
the  failure  envelope  and  the  point  ' P ' ,  defines  the  changing 
stability  of  the  system. 

By  fixing  the  angles  6  and  ^ t  the  plane  of  weakness  of 
the  material  is  determined.  The  variation  with  time  of  the 
distance  between  the  corresponding  point  P  on  the  surface  of 
the  Mohr's  circles  and  the  failure  envelope  will  result  in 
an  "INSTABILITY"  history  for  a  given  location  of  coordinates 
X,  Y,  Z.  This  risk  history  can  be  represented  as  a  curve  in 
an  instability  value  vs  time  diagram.  This  curve  is  then 
referred  to  as  an  instability  function. 
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Figure  7....  Definition  of  the  instability  function  in  three 
dimensions  as  the  distance  between  P  and  the  failure 
envelope . 
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This  kind  of  risk  function  can  be  used  to  study  induced 
seismicity  at  water  reservoirs  (Ur ibe-Carva jal  and  Nyland, 
1983),  in  such  a  way  that  this  function  will  depend  only  on 
the  loading  history  of  the  reservoir,  the  known  geological 
structures  (that  will  determine  the  angles  6  and  \ It),  and  the 
geometry  or  the  bathymetry  of  the  lake.  Stability  has  been 
defined  as  a  function  proportional  to  the  minimum  distance 
between  the  failure  envelope  and  'P'. 

The  use  of  a  Coulomb  failure  criterion  is  applicable 
when  rocks  behave  in  an  elastic  way  and  that  fracture  occurs 
in  a  brittle  way.  Although  rocks  behave  in  a  more 
complicated  way,  the  assumption  of  elastic  materials  is 
often  made  in  geophysics.  Solomon  et  al.  (1980)  and  many 
others  have  suggested  that  the  upper  few  tens  of  kilometers 
of  the  earth's  crust  can  be  treated  as  elastic  materials; 
Turcotte  (1974),  determined  that  the  upper  bound  for  this 
pseudo-elastic  behaviour  is  300°C.  This  temperature  acts  as 
the  limititation  to  the  applicability  of  this  kind  of  linear 
failure  envelope. 

Assuming  Mohr-Coulomb  failure  is  consistent  with  the 
assumption  that  the  incremental  stresses  cause  elastic 
deformation  particularly  near  failure,  the  assumption  of 
brittle  failure  may  not  be  true  for  all  faults,  but  it  is  a 
reasonable,  tractable  hypothesis. 

In  the  case  of  the  application  to  a  water  reservoir 
(chapter  5),  we  acknowledge  that  the  treatement  of  the 
as  porous  halfspace  consisting  of  an  elastic  matrix 
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saturated  with  water  is  a  simplistic  model.  However,  the 
stability  functions  are  relative,  and  only  serve  as 
indicators  of  how  the  risk  of  inducing  seismic  activity  is 
changing  with  respect  to  a  reference  initial  value. 

2.5  Residual  Stability  Function. 

The  determination  of  the  state  of  regional  stress 
itself  can  be  approached  by  determining  the  stress  field  in 
a  two  dimensional  finite  element  model.  A  constraint  on  this 
model  is  that  observed  seismicity  occurs  in  areas  of  high 
shear  stress.  These  areas  are  referred  to  as  being  unstable. 
The  stability  of  a  given  location  has  been  determined  by  the 
distance  of  a  Mohr  circle  from  a  Coulomb  type  of  failure 
envelope.  Since  the  values  obtained  represent  the  stress 
increment  that  is  required  to  reach  failure,  this  is  an 
instability  measure  (Ur ibe-Carva jal  and  Nyland,  1983). 

In  the  previous  section,  I  introduced  the  notion  that  a 
measure  of  the  seismic  stability  of  a  geologic  structure  was 
the  minimum  distance  between  the  envelope  describing  failure 
in  shear  and  the  Mohr  circle  for  the  stress  state.  This 
notion  is  applied  here.  A  large  value  for  this  quantity 
implies  a  low  risk  of  failure  or  high  stability.  Note  that 
stability  behaves  inversely  to  risk.  Another  geometric 
measure  of  the  relative  position  of  the  stress  with  respect 
to  the  failure  envelope  might  be  reasonable. 

I  introduce  in  this  part  a  measure  of  risk  of  failure  I 
call  residual  instability  (figure  8).  This  residual 
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Figure  8....  Shows  the  definition  of  residual  instability  as 
the  minimum  distance  between  the  Mohr  circle  and  the  Coulomb 
failure  envelope. 
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instability  is  derived  by  first  estimating  some  form  of 
average  stress  in  the  area  under  investigation  and  then 
determining  the  stress  in  an  anomalous  zone,  subtracting 
from  this  stress  the  ambient  background  stress  and 

\ 

calculating  the  distance  of  the  Mohr  circle  for  the  residual 
stress  from  the  failure  envelope.  The  normal  ( 1 i thostat ic ) 
stress  is  usually  that  which  would  result  from  a  gravitating 
uniform  half-space,  the  anomalous  stress  is  derived  from  the 
gravitational  stresses  that  would  be  present  when  a 
particular  geodynamic  structure  is  included.  In  chapter  3 
this  residual  instability  is  applied  to  the  study  of  seismic 
behaviour  and  the  determination  of  high  risk  areas. 

The  residual  instability  behaves  rather  differently 
from  our  previous  measure.  Residual  principal  stresses  can 
be  negative.  This  means  that  the  Mohr  circle  that 
corresponds  to  these  stresses  could  cross  a  failure 
criterion  that  intersects  the  origin  as  in  figure  8.  This 
does  not  mean  that  failure  will  occur,  it  merely  means  that 
the  distance  from  the  Mohr  circle  for  the  total  stress  to 
the  failure  envelope  has  decreased  and  failure  is  more 
likely.  Strictly,  failure  will  not  occur  until  the  Mohr 
circle  for  the  residual  stress  is  as  far  above  the  failure 
envelope  as  that  for  the  normal  stress  is  below  the  failure 
envelope.  The  residual  instability  can  be  negative.  This 
will  occur  when  the  Mohr  circle  for  the  residual  stress  is 
completely  below  the  failure  envelope.  What  is  important  is 
whether  this  residual  measure  is  larger  in  one  location  than 
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another  for  if  it  is,  the  location  with  the  larger  measure 
is  more  unstable. 

This  relative  treatment  of  the  problem  is  necessary. 

Although  the  assumption  that  the  strength  in  tension  of 

rocks  is  small  appears  justified  it  cannot  be  verified  on 

the  scales  we  consider  here.  The  sign  of  the  instability 

measure  depends  on  the  point  where  the  failure  envelope 

intersects  the  zero  normal  stress  axis  on  the  Mohr  diagram. 

Note  that  the  instability  values  considered  here  are 

distances  on  a  Mohr  diagram.  This  means  that  instability  has 

the  same  units  as  stress.  Keep  in  mind  however  that  an 

ambiguity  exists.  The  actual  intercept  of  the  failure 

envelope,  that  is,  the  strength  in  tension  of  the  formation 

is  an  unknown.  As  a  result  of  this  residual  instability 

values  for  any  model  are  nonunique  to  an  additive  constant. 

The  contours  of  equal  residual  instability  show 

specific  areas  where  the  resulting  stresses  combine  in  such 

a  way  that  it  is  possible  to  denote  them  as  minimum 
* 

stability  areas. 

2.6  The  Probability  Function. 

The  concept  of  instability  has  been  discussed 
previously  as  a  comparative  measure  between  two  states  of 
stress  and  represents  a  measure  of  the  risk  of  the 
occurrence  of  a  seismic  event.  It  deals  essentially  with  the 
determination  of  the  principal  stresses  and  its 
interpretation  in  form  of  a  Mohr  circle,  and  the  minimum 
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distance  between  this  circle  and  an  envelope  determined  by  a 
failure  criterion.  In  order  to  take  this  instability 
function  one  step  further  and  to  be  able  to  talk  in  a  clear 
language  about  risk,  I  introduce  here  a  change  to  the 
instability  or  risk  function  by  linking  it  to  the  probabilty 
of  failure  at  a  particular  location  of  a  model  with  respect 
to  the  rest  of  the  model  when  a  system  of  forces  is  applied 
to  it. 

This  probability  distribution  (figure  9)  is  also 
relative.  It  indicates  that  some  areas  will  be  closer  to 
failure  than  others  but  does  not  mean  that  at  a  high 
probability  valued  location  an  earthquake  must  occur.  It  is 
not  a  comparative  measure  of  stability  for  a  point  at  two 
different  times,  but  rather,  a  comparative  measure  for  all 
the  different  parts  of  the  model  at  a  particular  stage. 

To  calculate  the  probability  function,  first  it  is 
necessary  to  calculate  the  principal  stress  distribution  at 
every  point  of  the  model.  With  these  stresses  Mohr  circle 
diagrams  can  be  constructed.  In  figure  9  circle  1  is  the 
Mohr  diagram  that  corresponds  to  the  node  of  highest 
stability  or  minimum  risk,  which  as  seen  previously  is 
proportional  to  the  minimum  distance  between  the  Mohr  circle 
and  the  failure  envelope.  This  distance  can  be  negative  and 
in  the  figure  is  indicated  as  dm in.  Circle  2  corresponds  to 
that  node  for  which  the  resulting  distance  d  is  maximum. 

The  failure  criterion  specified  by  a  single  envelope  on 
a  Mohr-Coulomb  diagram  prescribes  an  absolute  condition  for 
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Figure  9 . . . . 


Definition  of  the  probability  function. 
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failure.  It  is  my  view  that  this  is  simplistic  for 
earthquakes.  I  prefer  to  think  of  a  probability  of  failure 
concentrated  in  the  neighbourhood  of  the  failure  envelope. 

In  what  follows  we  make  an  intuitive  argument  about  the  way 
the  finite  element  estimates  of  instability  could  be  mapped 
into  a  probability  function  describing  risk  of  failure. 

Stresses  can  be  calculated  at  every  node  in  a  finite 
element  model.  This  produces  a  set  of  numbers,  call  them 
{ d i }  for  the  distance  between  the  failure  envelope  and  the 
Mohr  Circle  at  the  ith  node  in  the  system.  Suppose  now  that 
each  value  d  can  be  mapped  onto  a  probability  distribution 
X(d)  in  a  unique  way. 

This  suggests  that  an  association  can  be  made  between 
the  position  of  the  Mohr  Circle  and  a  probability  of 
failure.  The  idea  is  illustrated  in  figure  9.  The  problem  to 
solve  is  the  determination  of  the  form  of  X(d). 

The  set  of  numbers  { d ; }  will  contain  both  small  and 
large  values.  It  need  not  be  obviously  a  set  of 
probabilities,  it  does  however  provide  a  link  to  the 
probabilities.  I  argue  first  that  a  normal  distribution  is  a 
reasonable  form  for  X.  All  probability  distributions  tend  to 
a  normal  distribution  near  their  peak  (the  central  limit 
theorem;  Reichl,  1980).  If  this  is  so  the  parameters  of  the 
distribution  must  be  determined.  In  particular  a  rule  must 
be  devised  whereby  a  probability  of  the  form 

y(z)  =  (27r)-l/2J^  exp  (-x  2  /a  2  )dx 


(2.10) 
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can  be  used  to  assess  the  risk  of  failure.  The  choice  of  the 
integrand  is  governed  by  the  requirements  of  a  normal 
distribution  and  the  integration  is  present  to  contain  the 
notion  that  the  risk  of  failure  is  cumulative  from  a 
threshold  z0.  Physically  one  can  think  of  the  failure 
criterion  as  being  controlled  by  the  density  of  cracks  of 
varying  size.  The  creation  of  a  catastrophic  cascade  crack 
is  described  by  a  probability  density  like  y(z). 

The  identification  is  complete  once  x,  o ,  z0,  and  z 
have  been  identified.  Let  dmx  be  the  largest  and  dm i n  the 
smallest  members  of  { d  f } .  Then  define  Rmx=dmx“dmin  and 
Ri=dmx“di .  Then  Xj=Ri/Rmx.  Since  the  variables  X;  range  from 
0  to  1  I  choose  o= 1  and  then  calculate  the  y  corresponding 
to  dj  from  the  integral 

y(z)  =  (27r)“l/2/z£  exp(  -x 2  )dx  (2.11) 

This  integral  is  easy  to  do  and  can  be  evaluated  rapidly 
when  Zi  is  between  0  and  1. 

Observe  that  the  integration  is  only  over  a  portion  of 
the  normal  distribution.  Choosing  the  lower  limit  of  the 
integral  at  z0  ensures  that  the  probability  of  a  seismic 
event  at  the  most  stable  node  is  zero.  This  is  to  some 
extend  arbitrary,  but  seems  justified  on  observational 
grounds.  Earthquakes  do  after  all  show  strong  clustering. 

The  choice  of  a  normal  distribution  is  now  dictated  by  the 
fact  that  we  are  'near'  the  high  probability  regions.  Since 
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we  never  integrate  over  the  entire  normal  distribution  no 
probability  will  become  1.  This  is  also  reasonable. 

A  question  now  arises  about  the  use  of  o= 1.  Since  the 
interpretation  done  is  relative,  I  compare  probabilities  at 
different  points,  the  value  of  o  does  not  really  matter.  An 
improved  approach  to  this  problem  could  involve  both  fitting 
o  to  data  from  models  and  examine  whether  other  than  normal 
distributions  could  work  better.  I  have  not  done  so  because 
normal  distributions  give  excellent  results,  as  will  be 
shown  later  in  this  thesis. 

Knopoff  (1971),  relates  the  probability  for  the 
occurrence  of  earthquakes  to  the  state  of  strain  of  the 
material.  However,  a  comparison  between  his  work  and  the 
definitions  presented  here  is  not  easy.  He  defines  a 
mathematical  algorithm  that  depends  on  three  different  types 
of  probabilities:  The  probability  that  the  stored  energy  of 
deformation  is  at  a  certain  level  P( E ft )dE ;  the  probability 
that,  if  this  energy  is  at  a  given  level,  an  earthquake  will 
occur  \(E);  and  the  transition  probability  that,  if  the 
earthquake  occurs  at  an  energy  state  X  the  final  energy 
state  will  be  at  a  given  level  E,  which  is  given  by 
T (X\E1dE .  These  are  used  to  predict  if  rupture  can  take 
place  at  a  particular  time  at  a  given  location. 

My  probability  function  could  be  seen  as  a  particular 
case  of  the  \(E)  of  Knopoff.  Here  I  assume  that  the  model  is 
at  a  given  energy  level.  Then  the  analysis  of  all  parts  of 
the  model  is  necessary,  which  is  a  static  exercise  highly 
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dependent  on  the  material  properties.  To  predict  changes  on 
the  distribution  of  probability  due  to  the  occurrence  of 
earthquakes  it  is  necessary  to  introduce  stress  drops  at  the 
focal  locations  and  recalculate  the  whole  model. 

It  should  be  noted  that  due  to  the  nature  of  this 
probability  distribution  a  zone  with  a  value  of  1  indicates 
where  the  model  is  most  unstable.  This  does  not  means  that 
an  earthquake  will  take  place  there,  but  rather  that  if  an 
event  should  occur  in  the  model  at  the  state  of  stresses 
modelled,  fracture  would  initiate  there. 

2.7  General  comments  on  the  instability  functions. 

Stability  functions  that  compare  all  the  parts  of  a 
model  at  a  particular  stage  (as  in  the  cases  of  the  residual 
instability  and  the  probability  functions),  can  be  used  to 
predict  the  probable  place  of  rupture  initiation.  To  be  able 
to  do  this  it  is  necessary  that  significant  differences 
exist  in  the  instability  distribution  within  the  model.  Once 
rupture  takes  place  the  crack  will  propagate  until  it 
reaches  a  zone  of  sufficient  stability  to  stop  this 
propagation.  This  will  occur  where  the  strength  of  the 
material  increases.  The  final  length  of  the  crack  could  then 
be  linked  to  the  length  of  the  high  risk  areas.  It  is 
possible  that  if  an  earthquake  occurs  it  might  propagate  to 
other  high  risk  areas  if  the  stable  bridges  between  these 
areas  and  the  one  in  which  the  event  took  place  are  narrow. 
These  kinds  of  predictions  are  demonstrated  in  chapters  3 
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and  4  of  this  thesis. 

In  instability  functions  dealing  with  several  stages 
will  be  possible  to  predict  in  which  stage  failure  might 
occur.  This  is  shown  in  chapter  5. 
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3,  Elastic  models  of  seismic  failure  and  its  application  to 

the  Valley  of  Mexico 

The  observation  that  near  boundaries  between  continents 
and  oceans  there  tend  to  exist  mountain  belts  (Gunn,  1947), 
that  in  general  are  parallel  to  subduction  arcs,  has  one 
exception  (Molnar  and  Sykes,  1969).  The  Mexican  volcano 
chain  known  as  the  Neovolcanic  axis  (which  is  the  mountain 
belt  associated  with  the  Middle  America  Trench,  north  of  the 
Tehuantepec  Ridge)  is  located  anomalously  far  back  from  the 
trench  and  is  not  situated  above  the  focus  of  intermediate 
depth  earthquakes.  The  strike  of  this  mountain  belt  differs 
by  about  15°  from  that  of  the  trench.  However,  it  is 
difficult  to  believe  that  the  volcanoes  are  unrelated  to  the 
subduction  zone.  Molnar  and  Sykes  (1969)  suggest  that 
factors  that  could  explain  this  anomalous  distribution  of 
volcanoes  may  be; 

1.  Magmas  in  this  region  are  not  generated  near  the  dipping 
seismic  belt. 

2.  Magmas  follow  an  indirect  path  to  the  surface  due  to  the 
existence  of  previous  zones  of  weakness. 

3.  The  slip  rate  of  the  downgoing  slab  may  have  changed 
during  the  last  few  million  years,  as  Mexico  approached 
the  East  Pacific  rise. 

This  third  hypothesis  might  lead  to  a  decrease  in 
seismic  activity  before  the  volcanism  stopped  (Molnar  and 
Sykes, 1969).  One  problem  with  this  mechanism  is  that  shallow 
activity  is  still  frequent  near  the  trench. 
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Other  factors  to  consider  in  search  of  the  probable 
mechanism  that  might  be  responsible  for  this  anomalous 
volcanic  distribution  is  the  absence  of  volcanism  between 
18°  and  15°  N  and  the  possibility  that  the  Tehuantepec  ridge 
may  have  been  active  in  the  past  10  million  years.  This  will 
be  dealt  with  in  chapter  6  of  this  work. 

In  this  chapter  the  residual  instability  function 
presented  in  section  2.5  is  applied  to  a  site  which  is 
seismically  active  and  could  be  related  to  the  Middle 
American  island  arc  north  of  the  Tehuantepec  ridge.  Thus  the 
behaviour  of  this  kind  of  instability  under  different 
circumstances  can  be  monitored  and  the  role  of  the 
subduction  mechanism  on  the  local  seismicity  detected  in 
thi's  location  can  be  tested.  This  chapter  is  an  expanded 
version  of  Elastic  models  of  seismic  failure  in  the  Valley 
of  Mexico  by  Ur ibe-Carva jal  and  Nyland  published  in  Pure  and 
Applied  Geophysics,  vol  120,  1982  in  which  it  was  proven 
that  the  regional  stresses  in  the  plate  have  little  effect 
on  the  distribution  of  shallow  unstable  zones  compared  to 
the  effects  of  local  geologic  structures  and  localized 
changes  of  stress. 

3 . 1  Introduction 

Restricted  parts  of  the  valley  of  Mexico  have  been 
plagued  by  swarms  of  small  earthquakes  during  recent  years. 
The  cause  of  these  events  is  puzzling,  but  they  may  relate 
to  massive  artificial  changes  in  the  hydrology  of  the 
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lake-bed  which  underlies  most  of  Mexico  City.  The  valley  is 
located  at  approximately  250  km  inland  of  the  Cocos-North 
American  plate  boundary. 

If  the  valley  responds  elastically,  plane  strain  finite 
element  calculations  can  be  used  to  relate  geologic  and 
seismic  structure  to  observed  seismicity  and  overall 
stability  of  the  valley.  These  calculations  show  that 
unstable  zones  and  areas  of  local  seismicity  coincide  if  the 
triggering  mechanism  of  the  seismicity  is  related  to 
variations  of  pore  pressure  and  water  load. 

In  order  to  make  such  a  connection  it  is  necessary  to 
ascribe  the  generating  mechanism  of  the  seismicity  to  the 
density  and  elastic  modulus  anomalies  associated  with  the 
valley  structure,  and  to  suggest  that  the  phenomenon  that 
triggers  this  seismicity  is  related  to  changes  in  water 
content  of  the  near  surface  formations.  Calculations  of  the 
change  in  pore  pressure  due  to  the  pumping  of  water  from  the 
system  show  that  the  unstable  areas  are  indeed  concentrated 
within  the  valley  and  in  general  tend  to  confirm  the 
conclusions  drawn  from  the  elastic  analysis. 

The  valley  of  Mexico  is  located  in  the  central  portion 
of  the  Neovolcanic  axis  of  Mexico  and  it  contains  Mexico 
City.  The  seismic  characteristics  of  the  valley  of  Mexico 
(Figueroa,  1971)  are  a  result  of  the  geological  conditions 
that  prevail  in  the  valley;  a  particularly  unfavourable 
geologic  condition  is  the  old  lake-bed  that  underlies  Mexico 
City  and  its  surrounding  areas.  The  seismic  behaviour  of  the 
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valley  is  characterized  by  high  shallow  microseismic 
activity  and  frequent  local  earthquakes  of  small  magnitude 
with  depth  in  the  2  to  8  km  range  (figure  10) 

This  seismicity  is  probably  a  symptom  both  of  the 
geological  processes  which  caused  the  valley  itself  and  the 
artificial  changes  in  the  area  since  the  beginning  of  the 
construction  of  the  city  (Oviedo,  1967).  In  the  valley  there 
has  been  and  still  is  volcanic  activity  and  there  are  also 
great  changes  in  hydrological  conditions.  These  processes 
have  occurred  continuously  and  have  created  multiple  faults, 
fractures,  volcanoes  and  barriers  that  have  closed  some  of 
the  natural  drainage  channels  of  the  valley.  As  a  result 
there  has  been  an  accumulation  of  a  variety  of  materials 
with  different  mechanical  properties  in  the  valley.  These 
heterogeneities  in  the  thickness  and  type  of  material  are 
responsible  for  the  variety  of  geological  processes  that 
have  been  observed.  Since  the  beginning  of  the  20th  century, 
the  expansion  of  the  city  and  the  lack  of  drainage  has  made 
it  necessary  to  pump  water  out  of  the  underlying  formations. 
This  is  a  probable  cause  of  the  observed  changes  of  the 
piezometric  levels  (Cruickskank  et  al.,  1979)  and  of  the 
sinking  suffered  in  some  areas  of  the  valley. 

Many  faults  have  been  located  from  geologic  evidence 
and  some  have  been  inferred  from  seismic  activity  but  many 
are  probably  covered  by  lava  flows  inside  the  valley  or  by 
the  layer  of  sediments  (which  in  some  locations  can  exceed  2 
km  of  thickness ) . 
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Figure  10....  Earthquake  locations  inside  the  Valley  of 
Mexico.  The  numbers  in  this  figure  indicate  known  locations; 
1  Ciudad  Universitar ia ,  2  Chapultepec , 3  Alameda  and  4 
Texcoco.  The  broken  line  outlines  the  valley  and  the 
continuous  one  Mexico  City 
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Any  investigation  of  the  geodynamics  of  the  valley  of 
Mexico  must  consider  the  probable  cause  of  the  local 
seismicity.  I  do  not  believe  that  it  is  directly 
attributable  to  volcanic  activity.  The  seismicity  is 
observed  in  the  South  Western  part  of  the  city  and  there  is 
no  evidence  of  thermal  activity  other  than  that  which  caused 
the  Chichinautz in  lava  flow  which  is  of  early  pleistocene 
age,  and  constitutes  the  southern  limit  of  the  valley.  There 
has  been  more  recent  activity  in  the  valley,  but  the  nearest 
recent  activity  is  about  70  km  East  of  the  area  in  which  the 
seismicity  was  observed.  The  region  of  Popocatepetl  is 
active  now.  A  connection  between  the  seismicity  observed  in 
South  West  Mexico  City  and  volcanic  manifestations  in  the 
South  East  is  probable  but  it  is  suggested  here  that 
triggering  by  pore  pressure  changes  due  to  pumping  of  water 
in  the  valley  is  a  more  plausible  explanation.  In  this 
chapter  the  relationship  between  the  geologic  and  hydrologic 

structure  of  the  valley  of  Mexico  and  its  seismic  behaviour 

* 

is  investigated. 

Little  is  known  about  the  geodynamic  characteristics  of 
the  valley.  The  importance  of  these  data  comes  from  the  fact 
that  the  seismic  activity  detected  inside  the  city  seems 
concentrated  at  a  few  locations,  and  that  the  damage 
suffered  by  the  structures  of  the  city  varies  greatly  in 
locations  only  a  few  hundred  meters  apart.  Methods  of 
obtaining  information  about  the  tectonic  state  of  stress  of 
the  valley  from  the  knowledge  of  its  geological 
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characteristics  will  be  investigated  in  this  chapter.  This 
stress  state  obviously  will  help  to  delineate  the  earthquake 
generating  mechanism. 

Solomon  et  al.  (1980)  pointed  out  that  the  state  of 
tectonic  stress  in  a  given  location  is  linked  to  the 
rheology  of  the  material  and  to  the  history  of  application 
of  a  system  of  forces.  None  of  this  information  is  known  in 
sufficient  detail  to  enable  the  determination  of  an  unique 
state  of  stress.  Nevertheless,  it  is  important  to  estimate 
this  stress  state  because  the  tectonic  stresses  cause 
crustal  deformation  and  control  many  tectonic  processes. 

In  this  study  an  elastic  model  of  the  crust  in  a  two 
dimensional  cross  section  will  be  used.  It  is  assumed  that 
the  only  force  acting  in  the  vertical  direction  is  the 
weight  of  the  material  itself.  This  elastic  treatment  is  in 
accordance  with  Lambeck  (1980),  who  pointed  out  that  the 
usual  elastic  theory  approximation  to  determine  stresses  in 
the  earth’s  crust  may  not  be  entirely  adequate,  unless  the 
region  in  study  is  sufficiently  far  away  from  the  plate 
boundaries  and  other  sources,  so  that  the  lithosphere  has 
cooled  and  elastic  forces  become  the  dominant  support 
mechanism  for  the  isostatic  load.  This  is  the  case  in  the 
valley  of  Mexico. 
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3„2  Previous  Sei smo-Geodynamical  Studies 

The  seismic  characteristics  of  the  valley  have  been 
recognized  since  the  begining  of  this  century.  Scientists 
have  focused  their  efforts  on  earthquake  problems  related  to 
seismic  engineering.  This  perhaps  is  the  main  problem  to 
overcome  when  more  general  geophysical  studies  are  intended. 
Although  the  area  has  been  monitored  continuously  since  at 
least  1909,  the  reports  are  more  concerned  with  intensities 
than  magnitudes,  the  location  and  reasons  for  the  events  are 
of  secondary  importance  to  the  effects  on  the  civil 
structures  of  the  city  and  the  determination  of  safe 
building  codes  for  the  different  areas  of  the  valley.  Only 
recently  have  geodynamical  studies  also  been  recognized  as 
also  of  importance.  This  is  understandable  if  one  considers 
the  human  factors  involved. 

All  this  implies  that  although  the  seismic  catalog  that 
is  available  is  complete  to  a  magnitude  threshold  of 
approximately  1.0,  the  magnitude  listing  in  the  catalog  is 
incomplete.  However,  the  original  smoked  paper  records  from 
the  Tacubaya  seismic  station  are  available  for  anyone 
interested  in  completing  the  task  and  from  these  at  least 
coda  magnitude  determinations  are  possible.  The  major 
problem  in  this  catalog  is  that  during  most  of  its  length 
the  only  seismic  station  available  is  Tacubaya,  which  makes 
the  epicentral  locations  and  the  depths  inaccurate. 
Nevertheless,  for  most  of  the  local  events  the  effects  of 
the  earthquake  in  the  valley  were  restricted  to  a  zone  of 
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small  radius  (1  km);  this  supports  the  locations  to  ±  2  km, 
and  the  epicentral  depths  for  these  local  events  can  also  be 
generalized;  they  are  shallow. 

The  seismic  instrumentation  of  the  Tacubaya  National 
Observatory  has  been  recording  without  interruption  since 
1909  (figure  11).  In  this  period  it  has  registered  more  than 
1900  local  events  inside  the  valley,  most  of  them  less  than 
25  km  from  the  observatory.  Figueroa  (1971)  calculated  that 
the  total  energy  released  in  these  local  earthquakes  is  of 
the  order  of  10' 7,5  ergs. 

Figueroa  (1971)  noted  that  a  common  feature  of  all 
these  events  is  a  dominant  short  period,  not  greater  than  .5 
sec,  independent  of  the  intensity  of  the  earthquake.  The 
time  of  duration  or  Coda  length  (T)  is  less  than  5  sec  for 
most  of  these  earthquakes.  This  results  in  a  magnitude  of 
0.5  at  the  most,  using  the  formula  (3.1)  for  coda  magnitude: 

M  =  A  -  B  log  T  (3.1) 

and  the  constants  A  and  B  determined  for  California 
(Richter,  1958).  A  magnitude  of  0.25  is  obtained  from 
applying  the  constants  for  Rangely  Colorado  (O'Neill  and 
Healy,  1973).  A  magnitude  of  0.4  is  the  result  when  using 
the  constants  determined  for  the  central  part  of  Chiapas  in 
southern  Mexico  (Ur ibe-Carva jal ,  1979).  These  calculations 
reinforce  the  previous  inference  that  the  Tacubaya  seismic 
catalog  is  complete  at  least  from  a  magnitude  of  1.0  for  the 
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Figure  11....  Number  of  events  recorded  in  the  Valley  of 
Mexico  vs  time  (modified  from  Figueroa,  1971). 
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whole  valley. 

Figure  11  shows  without  any  doubt  that  the  amount  of 
local  seismicity  recorded  inside  the  valley  has  been  slowly 
increasing  during  most  of  this  century.  It  is  reasonable 
then  to  assume  that  it  is  directly  related  to  the  growth  of 
Mexico  City  and  to  the  amount  of  water  pumped  out  of  the 
subsoil  in  the  area. 

The  region  has  been  monitored  for  more  than  70  years 
and  the  recurrence  times  for  large  earthquakes  in  the  Middle 
America  trench  zone  are  of  the  order  of  35  years.  If  the 
region  of  interest  is  a  tectonic  feature  of  the  subduction 
of  the  Cocos  plate  beneath  the  North  American  plate  it  can 
be  expected  that  the  recurrence  time  will  be  of  the  same 
order  in  the  valley.  However,  due  to  the  discrepancy  between 
the  strike  of  the  mountain  belt  and  the  trench,  the 
assumption  that  no  large  tectonic  earthquakes  occurred  in 
the  area  after  two  periods  is  weak.  Nevertheless,  no 
evidence  of  a  large  earthquake  is  found  in  the  past  200 
years . 

3.3  Seismic  instrumentation  inside  the  area 

As  mentioned  in  the  previous  section,  mos.t  of  what  is 
known  about  the  seismic  characteristics  of  the  valley  are 
due  to  the  records  obtained  at  the  Tacubaya  seismic 
observatory.  The  observatory  has  three  horizontal 
seismometers  and  one  Weichert  vertical  seismometer  with  a 
mass  of  17  tons,  a  period  of  1.5  sec,  and  a  maximum 
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amplification  of  2000  times.  All  this  intrumentat ion  belongs 
to  the  past  generation  of  mechanical  seismometers,  of  which 
there  are  few  still  in  operation.  Records  on  smoked  paper 
are  stored  and  any  event  of  interest  might  be  located  on 
them. 

In  1967  an  electromagnetic  seismic  station  began  to 
operate  in  the  old  Institute  of  Geophysics  of  the  University 
of  Mexico  (UNAM)  and  has  been  in  intermittent  operation 
since  then. 

Since  1961  the  Institute  of  Engineering  of  the  same 
university  has  set  a  network  of  seismic  equipment  throughout 
the  entire  valley.  By  1971  it  included  8  strong  motion 
recording  instruments  (accellerographs )  and  3  seismoscopes 
(figure  12)  (Figueroa,  1971). 

A  seismic  telemetering  network  was  installed  in  1976 
with  the  specific  goal  of  studying  the  seismic  behaviour  of 
the  valley.  In  cooperation  with  the  UNESCO  the  Institute  of 
Engineering  started  the  Sistema  de  Informacion 
Si smo-Telemetr ica  de  Mexico  (SISMEX).  By  1978  this  had  four 
seismic  stations  in  full  operation  (figure  12),  each  with  a 
short  period  (1  sec)  vertical  seismometer.  The  analog  signal 
from  these  stations  is  sent  to  a  central  location  in  the 
institute  itself  where  the  signals  are  recorded  on  magnetic 
tape  at  the  same  time  as  ink  galvanometer  seismograms  of 
each  station  are  drawn.  The  events  are  selected  from  these 
records  and  then  digitized  by  means  of  a  NOVA  800 
minicomputer.  The  data  are  stored  in  digitized  form  on 
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Figure  12....  Shows  the  location  of  the  seismic  stations  of 
the  SISMEX  network. 
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permanent  magnetic  tapes  which  allows  manipulation  for  any 
kind  of  seismic  study. 

This  network  has  been  integrated  into  a  nationwide 
network  since  1979,  which  allows  large  earthquakes  to  be 
studied  more  accurately.  The  local  events  in  the  valley  are 
rarely,  if  ever,  detected  by  seismic  stations  outside  the 
valley . 

The  SISMEX  network  has  improved  enormously  the  quality 
of  the  seismic  information  recorded  within  the  valley  and 
has  made  apparent  the  completeness  of  the  previous  catalog 
obtained  from  the  records  from  the  Tacubaya  seismic 
observatory  and  the  potential  for  high  accuracy  seismic 
studies  is  now  available.  This  network  has  operated  without 
interruption  and  information  derived  mainly  from  its  records 
is  now  being  analysed  to  understand  the  seismic  behaviour  of 
the  Valley  of  Mexico. 

3.4  State  of  stress 

The  determination  of  the  state  of  regional  stress 
itself  can  be  approached  by  determining  the  stress  field  in 
a  two  dimensional  finite  element  model.  A  constraint  on  this 
model  is  that  observed  seismicity  occurs  in  areas  of  high 
shear  stress.  These  areas  have  been  referred  to  as  unstable 
(Chapter  2).  The  stability  for  this  elastic  model  at  any 
given  location  can  be  determined  by  the  distance  of  a  Mohr 
circle  from  a  Coulomb  type  of  failure  envelope  as  shown  in 
figure  8  (Chapter  2).  The  precise  instability  criterion  used 
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in  this  chapter  will  be  the  Residual  Instability  Function  as 
described  in  section  2.5  of  this  thesis. 

To  be  able  to  calculate  the  desired  instability 
function  one  must  specify  a  model.  For  this  study  a 
seismically  determined  model  of  the  valley  of  Mexico 
(Havskov,  1982  b)  was  initially  adopted.  The  regional 
structure  of  the  area  is  modelled  as  a  set  of  horizontal 
layers  shown  in  figure  13.  This  model  was  obtained  from  an 
earthquake  swarm  that  occurred  in  February  1981,  and  is  an 
extension  of  one  based  on  refraction  profiles  (Havskov  and 
Singh,  1  978  )  . 

The  structure  of  the  upper  200  meters  of  sediments  is 
well  known  since  numerous  wells  have  been  drilled  in  the 
valley.  Unfortunately,  only  one  has  gone  to  a  depth  of  2000 
meters  (Oviedo,  1967).  However,  since  the  regional  state  of 
stress  is  the  main  concern  here,  the  details  of  this 
structure  are  of  secondary  importance.  The  general 
conclusions  of  detailed  studies  (Cruickshank  et  al.,  1979) 
will  be  taken  to  fix  the  average  behaviour  at  shallow 
depths . 

The  models  will  deal  at  the  most  with  the  upper  20  km 
of  the  lithosphere.  For  the  elastic  models  a  finite  element 
program  that  works  with  linear  elements  of  three  and  four 
sides  has  been  developed.  The  results  obtained  have  been 
reproduced  also  by  applying  the  program  ADINA  ( Bathe , 1 978 ) . 

A  general  description  of  ADINA  can  be  found  in  chapter  4. 

The  porous  media  analysis  was  done  with  the  help  of  the 
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program  ADINAT  (Bathe,  1981)  which  is  a  complementary  finite 
element  program  to  ADINA  allowing  the  modelling  of 
consolidation  problems. 

The  grids  of  finite  elements  used  in  the  initial  models 
have  in  common  that  for  the  upper  6  km  the  elements  taken 
are  cubes  of  approximately  1  km  per  side  while  below  this 
depth  the  elements  have  a  dimension  of  2  km  in  the  vertical 
direction.  This  rule  was  changed  for  the  final  models  due  to 
the  complicated  geology  and  the  node  distribution  used  can 
be  seen  in  figure  16a. 

The  initial  model  (figure  13)  is  a  modification  of 
Havskov's  model  in  the  sense  that  the  horizontal  layers  have 
been  placed  inside  half  of  a  valley  shaped  structure.  The 
dip  of  the  walls  in  the  valley  was  initially  taken  as  45°. 
The  variation  of  the  stability  function  can  then  be  observed 
as  a  function  of  this  angle. 

The  boundary  conditions  taken  are  as  follows:  on  the 
bottom,  no  vertical  displacement  can  take  place.  The  valley 
is  taken  to  be  symmetrical,  that  is,  the  model  is  restricted 
to  move  in  the  vertical  direction  only  on  the  right  side.  On 
the  left,  a  horizontal  stress  of  0.8  of  the  gravitational 
load  is  applied.  The  value  0.8  is  used  because  the  maximum 
shear  stress  calculated  with  this  criterion  gives  values 
similar  to  those  determined  by  the  maximum  shear  stress 
curve  for  hard  rock  (McGarr,  1980),  and  those  determined 
from  hydrofracturing  experiments  by  Haimson  (1978).  It  is  in 
agreement  with  the  few  available  measurements,  from  which 
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Figure  13....  a)  Shows  the  resulting  instability 
distribution  in  bars,  b)  Shows  at  the  bottom,  the 
displacement  suffered  due  to  the  gravitational  load. 
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McGarr  (1980)  indicated  that  the  shear  stress  increases  with 
depth  in  the  upper  3  to  5  km,  and  there  is  no  sign  of  a 
reduction  in  its  gradient  with  depth. 

The  elastic  moduli  used  were  calculated  by  fixing  a 
Poisson  ratio  of  .25  for  all  the  layers  and  giving  the 
following  densities  for  each  of  the  strata  from  above 
downwards:  2.2,  2.6,  2.65,  and  2.8  g/cm3  for  the  basin 
structure . 

Figure  13c  shows  the  relation  between  the  different 
components  of  the  stress  tensor  and  the  resulting  residual 
instability  distribution  for  a  model  similar  to  the  previous 
one.  The  same  densities  are  used  in  this  exercise,  however, 
the  boundary  conditions  at  the  vertical  walls  of  the  model 
differ  from  the  other  models  used  in  this  chapter.  The  nodes 
at  these  boundaries  can  not  move  in  the  horizontal 
direction.  That  is,  the  only  force  that  is  acting  in  the 
model  is  gravity. 

The  measure  defined  as  residual  instability  (chapter  2) 
will  be  explored  in  this  chapter.  The  residual  stress  will 
be  derived  from  the  change  in  gravitational  stresses  that 
results  from  the  introduction  of  the  structure  of  the  Valley 
of  Mexico  in  the  homogenous  halfspace.  Detailed  observations 
of  the  regional  stress  in  the  Valley  of  Mexico  are  not 
available  but  it  is  known  that  thrust  earthquakes  have  been 
observed  there  (Havskov,  1982. b).  The  present  model  permits 
such  events  because  of  the  inclusion  of  horizontal 
compression  by  constraining  the  boundaries  of  the  finite 
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Figure  13c....  Stress  components  and  the  corresponding 
residual  instability  distribution  for  a  gravitational  model 
of  the  Valley  of  Mexico. 
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element  model. ^ 

The  nature  of  faulting  in  this  kind  of  model  depends  on 
the  relative  magnitudes  of  the  principal  stresses.  Generally 
in  tectonics  one  is  vertical,  denoted  a3.  If  a3  is  greater 
than  both  horizontal  stresses  a,  and  a2r  the  regime  is  in  a 
normal  faulting  condition.  This  is  the  case  for  the  region 
of  study  if  the  anomalous  structure  of  the  valley  is  not 
included.  If,  however,  a  low  density  inclusion  such  as  a 
valley  is  placed  in  the  model  the  stress  variation  induced 
by  gravitation  can  lead  to  thrust  faulting.  This  is 
particularly  true  if  in  addition  to  the  gravity  the  stress 
field  is  modified  by  external  agencies  such  as  regional 
stress  fields  which  can  be  modelled  as  boundary  conditions 
on  a  finite  element  model. 

The  contours  of  equal  residual  instability  in  figure 
13a  show  minimum  residual  stability  areas.  The  displacement 
of  each  of  the  nodes  of  the  model  is  given  by  the  distortion 
of  the  grid  of  figure  13b.  A  compression  of  the  material  in 
the  zone  of  low  stability  has  occurred  due  to  the 
concentration  of  stress  resulting  from  the  density  contrast 
between  materials.  These  displacements  have  been  increased 
by  a  factor  of  5  in  order  to  make  the  change  visible. 

By  comparing  figure  13a  with  figure  14  it  can  be  seen 
that  as  the  slope  of  the  valley  wall  increases,  the  zone 
located  beneath  the  valley  becomes  more  stable.  The 
sediments  of  the  valley  behave  in  a  more  complex  fashion  but 
usually  show  greater  instability  for  steep  sided  valleys 


1 


85 


than  for  shallow  slopes  on  the  valley  walls  (figure  14). 

To  incorporate  an  external  stress  in  a  regional  sense, 
it  should  be  recalled  that  the  driving  tectonic  stress 
introduced  by  the  movement  of  the  plates  is  in  the  range  of 
200  to  300  bars  (Solomon  et  al,  1980).  This  can  be  modelled 
by  increasing  the  stress  by  400  bars  in  the  horizontal 
direction  applied  at  the  left  side  of  the  models.  No 
significant  change  was  observed  in  the  location  of  the 
unstable  zones.  Therefore,  instability  is  mostly  dependent 
on  the  density  contrast,  that  is,  on  the  geometry  of  the 
model  itself. 

It  is  generally  thought  that  there  might  be  one  or 
several  volcanic  cones  buried  beneath  the  sediments.  For 
this  reason  a  model  containing  a  cone  in  the  bottom  of  the 
valley  was  studied.  The  resulting  instability  distribution 
shown  in  figure  15a,  resembles  that  of  figure  13a  with  some 
minor  changes:  the  deeper  parts  of  the  valley  are  slightly 
more  stable,  the  zone  beneath  the  slope  has  stabilized 
somewhat  but  the  zone  beneath  the  valley  has  become  less 
stable.  In  order  to  approach  the  problem  in  a  more  realistic 
way,  the  smooth  wall  of  the  valley  was  roughened;  its 
stability  distribution  is  shown  in  figure  15b.  The  unstable 
zones  from  figures  13a  and  15a  are  somewhat  enhanced  (figure 
15)  . 

Several  things  can  be  learned  from  these  simple  models. 
One  is  that  by  introducing  structures  it  is  possible  to 
change  the  shape  of  the  instability  curves  in  such  a  way 
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Figure  14....  Shows  the  residual  instability  distributions 
in  bars,  resulting  from  varying  the  slope  angle,  a)  For  a 
slope  of  60°  and  b)  For  a  slope  of  30°. 


87 


i _ _ _ » 

K  IV1 


Figure  15....  Shows  the  residual  instability  distribution  in 
bars  for  the  following  variations:  a)  a  cone  has  been 
introduced  in  the  bottom  b)  the  left  side  of  the  valley  has 
roughened . 
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that  the  maximum  values  correspond  to  the  areas  where 
earthquakes  have  been  detected.  This  brings  in  the  danger  of. 
exaggerating  this  procedure  in  order  to  force  a  result  in  a 
poor  initial  model.  Clearly,  to  reduce  the  inherent 
nonuniqueness  it  is  necessary  to  support  the  conclusions 
with  other  information. 

In  order  to  improve  the  model  in  the  upper  3  km  of  the 
crust,  use  was  made  of  the  geologic  map  of  the  valley 
(Instituto  de  Geologia,  1968).  Beneath  that  depth,  the 
Havskov  model  was  adopted  for  the  final  model  analysed 
(figure  16).  For  this  analysis  the  same  boundary  conditions 
were  used  as  in  the  previous  models. 

Figure  16b  shows  a  stable  zone  bounded  by  contours  of 
value  .15  at  a  depth  of  3  to  5  km.  This  stable  zone  appears 
to  be  discontinuous  in  the  stippled  area  which  is  a  section 
through  the  seismically  active  volume.  The  stippled  area 
covers  the  depths  at  which  events  occur  and  is  the  part  of 
the  seismic  zone  indicated  in  the  plan  on  figures  17  and  19 
which  lies  under  the  line  AB.  Figure  16b  also  shows  an 
instability  area  below  6  km  which  suggests  a  potential  for 
seismic  events  at  depths  below  the  aquifer. 

The  results  of  the  preliminary  models  shown  in  figures 
13  to  16  all  show  that  the  zone  at  the  bottom  of  the 
sediments  at  the  flank  of  the  basin  (lower  corner  of  the 
valley)  introduces  a  concentration  of  stresses  that  results 
in  low  stability  areas.  By  considering  the  valley  as  a 
region  inside  a  circle  of  mountains  a  low  stability  ring 
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Figure  16....  Shows  a)  the  nodes  used,  and  b)the  residual 
instability  distribution  for  a  cross-section  on  the  line  AB 
(figure  17).  The  shaded  area  represents  the  location  of 
observed  seismicity. 
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around  the  valley  is  delineated,  the  ring  appearing  at 
different  positions  at  different  depths.  However,  estimates 
of  depths  of  seismicity  are  difficult  to  obtain  with  the 
best  of  seismic  networks,  and  since  the  data  discussed  here 
were  obtained  with  less  than  ideal  observing  networks  it  is 
considered  that  the  data  are  not  accurate  enough  to  reliably 
discriminate  mechanisms  with  depth. 

Other  than  the  data  for  the  1981  (Havskov,  1982. b) 
swarm  which  is  located  to  ±1  km,  supporting  data  (Figueroa, 
1971)  rely  on  historical  evidence  and  damage  reports.  These 
are  sufficient  to  isolate  the  events  to  the  region  indicated 
on  figures  17  and  19  but  not  sufficient  to  publish  more 
precise  locations.  Nevertheless  the  triggering  mechanism 
probably  does  vary  with  depth  and  some  of  the 
inconsistencies  between  the  hydrologic  behaviour  and  the 
model  could  be  explained  this  way.  The  ring  defined  by  the 
depth  of  the  bottom  of  the  valley  is  considered  to  be 
representative.  The  area  where  earthquakes  have  been 
observed  (figure  17)  falls  entirely  within  this  ring. 

3.5  Discussion  of  Elastic  Models 

In  order  to  link  seismicity  and  the  instability  zones 
indicated  by  the  finite  element  model  calculations  it  is 
necessary  to  propose  a  mechanism  that  might  trigger  events 
in  the  anomalous  zone  of  figure  16b  only  at  the  location  of 
the  observed  seismicity.  Why  events  occurred  only  in  a 
restricted  part  of  the  stability  low  must  be  explained.  Here 
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it  is  suggested  that  the  triggering  mechanism  could  be 
related  to  pumping  induced  subsidence. 

Cruickshank  et  al.(1979)  developed  a  finite  element 
model  of  the  upper  crust  of  the  valley  and  were  able  to 
determine  the  thickness  of  the  water  saturated  layer,  the 
thickness  of  the  layer  over  it,  the  change  of  piezometric 
levels  and  the  rate  of  sinking  of  the  valley.  They  stated 
that  there  is  a  direct  relationship  between  the  observed 
sinking  and  the  pumping  of  water  from  the  water-saturated 
layer  and  that  the  only  way  of  stopping  the  subsidence  would 
be  to  stop  the  pumping.  Their  model  of  the  valley  consisted 
of  a  water  saturated  aquifer  overlying  the  basement  and 
underlying  a  narrow  non-saturated  layer.  This  model  is  an 
approximation  to  the  detailed  model  studied  by  Marsal  and 
Mazari  (1969).  Here,  the  upper  layer  was  treated  as  perched 
aquifers  separated  by  thin  impermeable  strata.  Marsal  and 
Mazari  pointed  out  that  as  a  consequence  of  the  exploitation 
of  ground  water,  a  strong  alteration  of  the  piezometric 
levels  had  occurred.  From  the  nature  of  equi-pressure  drop 
contours  in  the  area  (figures  XII-10  and  XI 1—11  from  Marsal 
and  Mazari,  1969),  they  inferred  that  the  effect  in  the 
aquifers  was  more  pronounced  in  the  west  section  of  the 
city.  Here  they  found  that  not  only  the  gradients  are 
larger,  but  also  maximum  pressure  drop  occurs.  Using  these 
facts  they  proposed  that  the  piezometric  level  drops  are 
related  to  the  underground  recharge  and  concentration  of 
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The  results  of  Cruickshank  et  al.  (1979)  predict  the 
sinking  and  the  piezometric  change  in  the  valley  accurately. 
Therefore  for  the  regional  analysis  proposed  here  and  due  to 
the  reasons  explained  earlier,  the  model  proposed  by 
Cruickshank  et  al.  will  be  used  in  the  search  for  a 
triggering  mechanism.  If  in  fact  several  aquifers  exist, 
they  must  be  connected  by  many  multiple  faults  as  suspected 
by  Figueroa  (1971).  The  rate  of  change  of  the  piezometric 
drop  in  the  upper  formations  was  found  to  have  a  high 
between  the  observed  values  of  20  and  30  (Cruickshank  et 
al.,  1979)  (figure  17).  This  value  represents  the 
piezometric  drop,  measured  as  the  difference  of  readings 
between  a  piezometer  placed  in  the  aquifer  and  one  placed  in 
the  water  table,  and  is  given  in  meters.  The  location  of  the 
earthquakes  mentioned  above  coincides  with  the  area  where 
the  gradient  of  this  piezometric  drop  is  high.  In  this  area 
the  velocity  of  flow  may  be  expected  to  increase. 

The  zone  where  the  pressure  in  the  pores  of  the  rock 
has  changed  the  most  is  where  the  effective  stress  has 
suffered  the  greatest  change.  The  effective  stress  that 
describes  the  behaviour  of  porous  media  is  the  result  of  the 
difference  between  the  normal  stress  and  the  pore  pressure. 
It  has  been  shown  by  several  authors  e.g.  Gough  (1978), 
Simpson  and  Negmatullaev  (1981),  Withers  and  Nyland 
(1976,1978)  and  others,  that  the  introduction  of  a  water 
load  can  change  the  state  of  stress  in  a  given  area  and  can 
induce  some  seismic  activity  depending  on  the  geological 
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Figure  17....  Observed  piezometric  drop  in  the  valley  given 
in  meters.  The  shaded  area  represents  the  location  of 
observed  seismicity. 
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conditions  prevailing  in  the  area.  In  the  case  of  artificial 
lakes,  the  induced  activity  seems  to  be  concentrated  in  the 
upper  10  km  of  the  earth's  crust. 

Although  the  general  belief  is  that  the  increasing 
water  load  will  decrease  the  effective  stress  in  such  a  way 
that  induced  seismicity  can  occur,  decreasing  a  load  may 
also  trigger  events  in  an  area  close  to  failure  (figure  18). 
This  initial  effect  can  be  reduced  by  diffusion  of  pore 
pressure  if  the  area  is  not  in  a  thrust  faulting 
environment.  If  an  increase  in  water  load  will  increase  the 
risk  of  seismic  activity  with  normal  faulting  and  tends  to 
stabilize  an  area  prone  to  thrust  faulting  (Gough,  1978  and 
Withers  and  Nyland,  1978),  the  decrease  of  water  load  will 
act  oppositely  in  the  different  stress  environments.  Figure 
18  demonstrates  that  the  reduction  of  water  load  IN  the 
aquifer  will  in  fact  reduce  the  value  of  the  stability 
function  in  the  zone  located  BENEATH  the  aquifer.  This  is  so 
because  in  this  region  the  major  effect  of  pumping  water  is 
a  reduction  of  the  normal  stress  on  a  less  permeable  stratum 
below  the  aquifer.  Since  the  effect  is  due  to  load  changes 
it  can  be  expected  to  occur  rapidly. 

Such  changes  have  two  effects  on  the  state  of  stress  on 
the  region;  one  is  the  increment  of  the  stress  due  to 
changes  in  the  mass  of  the  aquifer,  the  other  is  the  change 
of  the  effective  stress  due  to  a  delayed  change  of  pore 
pressure  at  depth  (Gough,  1978).  These  phenomena  can  be 
studied  with  Mohr  diagrams  which  can  be  plotted  for  normal 
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Figure  18....  Stability  changes  due  to  the  change  in  the 
gravitational  load  for  different  stress  environments.  a) 
normal  faulting,  b)  thrust  faulting  and  c)  strike-slip 
faulting  enviroment. 
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faulting  regimes  (figure  18a),  thrust  faulting  regimes 
(figure  18b),  and  strike  slip  faulting  regimes  (figure  18c). 
It  should  be  recognized  that  in  each  of  these  regimes  the 
maximum  and  minimum  principal  stresses  are  oriented 
differently  and  as  a  consequence  the  application  of  a  fluid 
load  will  cause  different  increments  in  the  Mohr  circle. 

Circles  1  in  figure  18  represent  the  initial  state  of 
stress  while  circles  2  represent  stresses  resulting  from 
subtracting  the  weight  of  the  water  above  the  stress 
distribution  environments.  It  is  clear  from  these  diagrams 
that  unloading  decreases  stability  immediately.  Although 
pore  pressure  diffusion  will  tend  to  stabilize  the  result 
the  impermeable  base  of  the  aquifer  may  isolate  the  unstable 
zone  from  this  phenomenon. 

It  is  likely  that  the  actual  sinking  is  representative 
of  the  amount  of  water  removed  at  the  different  locations 
inside  the  valley,  so  lines  of  equal  sinking  have  been 

overlaid  on  the  area  of  occurrence  of  the  earthquakes 

✓ 

(figure  19).  Clearly,  subsidence  alone  is  insufficient  to 
explain  the  presence  of  events  only  at  one  part  of  the  low. 

The  fact  is  that  pumping  of  water  has  been  stronger  in 
recent  years  in  the  low  stability  area  (Cruickshank , 
personal  communication),  which  means  that  the  overall 
sinking  reported  is  not  completely  consistent  with  the 
amount  of  water  retrieved  at  particular  locations.  The 
suggestion  that  the  amount  of  water  pumped  out  in  this  zone 
may  be  a  triggering  mechanism  is  supported  by  the  thrust 
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Figure  19....  Shows  the  observed  subsidence 
given  in  meters.  The  shaded  area  represents 
observed  seismicity. 


of  the  valley 
the  location  of 


98 


mechanisms  determined  for  the  1981  earthquake  swarm  by 
Havskov  (1982  b) .  This  is  so  because  a  reduction  in  pore 
pressure  can  activate  a  region  close  to  thrust  failure. 

The  results  obtained  by  Cruickshank  et  al. (  1979) 

% 

suggest  that  the  stability  is  decreased  far  more  by  the 
change  in  the  effective  stress  in  the  area  where  the  zones 
of  low  stability  and  the  high  rate  of  piezometric  drop 
intersect  (figures  17  and  19).  The  high  velocity  of  the  flow 
and  the  high  rate  of  change  of  the  amount  of  water  in  the 
overlying  formations  affect  the  effective  stresses  by 
introducing  load  changes  in  the  underlying  formations  and 
pore  pressure  fluctuation  in  the  aquifer. 


3.6  Porous  media  modelling 

It  has  been  recognized  for  a  long  time  that  the 
increase  of  water  content  in  a  geological  environment  can 
trigger  earthquakes.  This  is  possible  if  the  area  is  prone 
to  normal  or  stike-slip  faulting  mechanisms.  This  increase 
in  water  content  could  be  due  to  the  impounding  of  an 
artificial  lake  or  to  the  direct  injection  of  fluids  by 
boreholes . 

The  evidence  for  the  second  type  of  induced  seismicity 
comes  mostly  from  the  Denver  and  the  Rangely  fluid  injection 
studies.  Healy  et.  al.  (1968)  were  able  to  detect  in  the 
Denver  experiment  a  close  relationship  between  the  times  of 
injection  of  waste  fluids  and  the  occurrence  of  seismic 
activity  in  the  surrounding  area.  They  were  able  to  infer 
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that  the  activity  is  delayed  approximately  one  year  from 
when  the  injection  took  place.  The  epicenters  determined  lay 
within  an  elliptic  zone  of  10  km  by  10  km  which  encloses  the 
disposal  well.  The  events  recorded  have  a  strike-slip 
mechanism  which  is  consistent  with  the  faults  observed  in 
the  area  and  have  magnitude  up  to  5.  There  is  no  evidence  of 
previous  seismic  activity  in  the  area.  The  well  was  drilled 
through  3700  m  of  sedimentary  rocks  and  finished  in 
fractured  precambrian  granitic  gneiss. 

At  the  Rangely,  Colorado  Oilfield  the  U.  S.  Geological 
Survey  carried  out  an  experiment  to  control  the  seismic 
activity  of  the  area  (Raleigh  et  al,  1976).  They  were  able 
to  verify  the  effects  of  pore  pressure  changes  on  the 
effective  stresses,  and  relate  these  stresses  determined 
from  field  observations  to  the  occurrence  of  seismic  events 
using  a  Mohr-Coulomb  failure  criteria.  The  oilbearing 
formation  is  Weber  sandstone  which  is  350  m  thick  and  at  a 
depth  of  1700  m.  This  overlays  other  sandstone  and  limestone 
formations  of  Paleozoic  age.  The  basement  is  found  at  a 
depth  of  3000  m.  The  induced  events  had  a  strike-slip 
mechanism  which  is  also  in  accordance  with  the  observed 
fault  pattern  of  the  area. 

In  the  previous  section  it  has  been  suggested  that  if 
delayed  earthquakes  with  normal  mechanisms  can  be  induced  by 
increasing  the  volume  of  water  involved,  it  is  also  possible 
to  induce  instantaneous  seismic  events  with  thrust 
mechanisms  by  decreasing  it.  This  was  done  by  splitting  the 
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response  of  the  media  into  two  parts,  a  rapid  reaction  due 
to  the  weight  of  the  fluid  alone  and  a  time  delayed  effect 
due  to  pore  pressure  relaxation. 

Since  the  models  presented  to  this  point  did  not 
consider  pore  pressure  changes,  the  instability  values  could 
be  overestimated  by  these  calculations.  To  test  this 
hypothesis  a  model  consisting  from  porous  materials  was 
studied  (figure  20).  Several  assumptions  were  made: 

1.  The  different  materials  described  in  all  the  other 
models  are  now  taken  to  be  porous  materials,  consisting 
of  elastic  matrices  with  the  voids  saturated  with  water. 

2.  Permeabilities  were  assigned  to  each  of  these  materials 
and  were  chosen  in  such  a  way  that  they  were 
representative  of  the  type  of  material  and,  as  far  as  is 
known,  of  their  geologic  characteristics.  The  material 
properties  assumed  in  our  model  are  in  accordance  with 
the  geologic  data  of  the  area,  and  the  vertical 
distribution  of  seismic  velocities  observed  there. 

3.  In  order  to  allow  for  a  finite  element  solution  the 
boundary  conditions  were  taken  in  such  a  way  that  no 
water  can  go  through  the  lower  boundary  of  the  model. 
This  means  that  all  the  water  is  contained  within  the 
upper  part  of  the  earth's  crust;  in  our  model  the 
impermeable  stratum  is  found  at  12  km  of  depth.  This  was 
chosen  since  the  main  area  of  interest  here  is  the 
behaviour  of  the  top  layers  of  sediments. 
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Figure  20....  Shows  the  change  in  pore 
bars,  calculated  due  to  the  pumping  of 
underlying  formations  of  the  Valley  of 


pressure  given 
water  from  the 
Mexico . 
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The  permeabilities  assigned  are:  30.0,  0.5,  4.0,  2.0, 
and  0.9  millidarcys  for  the  layers  taken  from  top  to  bottom 
of  figure  20.  The  extraction  of  water  is  modelled  by  placing 
suitable,  designed  extractors  at  the  surface  nodes.  The  rate 
at  which  these  "fluid  pressure  sinks"  remove  water  was 
adjusted  to  reproduce  the  observed  subsidence  in  the  Valley 
of  Mexico.  There  were  two  phases  of  constant  rate  of 
pumping.  The  first  60  years  accounted  for  90%  of  the  total 
load  in  the  part  located  to  the  right  of  the  high  risk  area 
of  figure  18b  and  no  pumping  was  done  in  the  high  risk  area. 
During  the  second  phase,  lasting  5  years,  the  pumping  rate 
was  such  that  the  total  load  applied  during  the  65  years  was 
equivalent  in  volume  to  the  observed  subsidence  at  each 
location . 

The  isolines  of  pore  pressure  drop  that  result  from 
this  analysis  should  be  interpreted  as  the  long  term 
potential  for  relaxation  of  the  stresses.  They  indicate  how 
much  the  effective  stress  at  a  particular  location  can 
change  due  to  future  pore  pressure  relaxation.  Greater 
values  on  figure  20  suggest  that  the  location  will  become 
more  stable  once  the  pore  pressure  effect  is  introduced. 

To  be  able  to  incorporate  this  porous  media  analysis, 
to  the  overall  instability  distributions  presented  in  this 
chapter,  it  must  be  noted  that  figures  13  to  16  represent 
the  state  of  stress  after  the  impounding  of  the  water  alone. 
They  do  not  consider  any  change  of  stress  due  to  pore 
pressure  variation.  With  this  in  mind,  the  combination  of 
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the  results  presented  in  figures  16  and  20,  will  represent  a 
measure  of  long  term  instability.  The  combination  of  a  high 
residual  instability  value  (figure  16b),  with  an  area  of 
large  (positive)  change  of  pore  pressure  (figure  20) 
indicates  that  an  unstable  area  will  not  recover  its 
stability  until  after  a  certain  time  period  unless  external 
stresses  are  introduced.  This  shows  the  potential  of  these 
areas  to  be  high  risk  value  zones.  One  of  the  possible 
external  sources  of  stresses  are  earthquakes. 

Comparing  figures  16b  and  20,  it  becomes  apparent  that 
the  shallow  seismicity  will  diminish,  if  not  disappear, 
after  the  effect  of  pore  pressure  changes  is  assimilated. 
This  strongly  suggests  that  the  events  are  directly  related 
to  the  extraction  of  the  water  and  will  disappear  some  time 
after  the  pumping  stops.  However,  the  unstable  areas  below  2 
km  of  depth  coincide  with  the  zones  at  which  the  pore 
pressure  attains  greater  values  at  the  same  depths.  This  can 
be  interpreted  to  mean  that  if  an  earthquake  occurs  at  these 
depths  it  will  be  located  in  these  areas. 

It  is  understood  that  the  treatment  of  the  lithosphere 
as  a  porous  material  consisting  of  an  elastic  matrix 
saturated  with  water  is  a  highly  idealized  model  but  it  is 
used  here  to  test  if  in  fact  the  seismicity  of  the  area 
could  be  linked  to  the  change  in  hydrologic  conditions.  The 
use  of  more  complicated  models  will  not  affect  this  result. 

It  should  also  be  pointed  out  that  although  this  study 
and  the  Denver  and  Rangely  experiments  are  related  to 
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induced  seismicity  several  differences  between  them  exist. 

At  the  Valley  of  Mexico  the  earthquakes  observed  have  thrust 
faulting  mechanisms,  the  activity  seems  to  have  increased  in 

relation  to  the  higher  rate  of  water  EXTRACTION  while  Denver 

\ 

and  Rangely  had  seismic  events  with  strike-slip  mechanisms 
delayed  with  respect  to  the  time  of  fluid  INJECTION.  In  all 
three  cases  the  earthquake  mechanisms  are  in  accordance  with 
the  faulting  pattern  observed  in  each  area. 


3.7  Conclusions 

The  models  used  for  the  stability  determinations  are 
restricted  to  elastic  materials  and  do  not  take  into 
account,  except  in  a  qualitative  way,  diffusion  of  water 
through  the  media.  As  a  result,  the  curves  in  the  upper 
layer  of  figure  16b  are  only  approximate  (an  alluvial 
formation).  Beneath  this  depth  the  results  are  credible.  The 
general  exercise  is  model  dependent  however,  and  the  model 
is  far  from  perfect.  Nevertheless,  the  general  features  of 
the  stability  distribution  should  be  independent  of  small 
changes  in  the  model. 

These  studies  have  shown  that  from  qualitative  geologic 
knowledge,  and  some  knowledge  of  induced  seismicity  due  to 
the  change  in  pore  pressure,  something  can  be  learned  about 
the  nature  of  the  seismicity  observed  in  the  valley  of 
Mexico. 

The  following  conclusions  have  been  reached: 

1.  The  concentration  of  stresses  due  to  the  bottom  corner 
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of  the  valley  creates  a  low  stability  zone  which 
includes  the  zone  where  seismic  activity  has  been 
observed. 

2.  There  seems  to  be  a  relation  between  the  observed  events 
and  some  parts  of  the  minimum  stability  zone  described 
in  figure  16b. 

3.  The  triggering  mechanism  in  those  parts  of  the  minimum 
where  the  activity  is  located  might  very  well  be  related 
to  the  hydrological  conditions  of  the  valley. 

4.  The  triggering  mechanism  of  the  deeper  events  may  be 
related  to  the  decrease  in  water  content  in  the  aquifer. 
There  may  also  be  other  tectonic  processes  at  these 
depths,  or  the  aquifer  may  not  be  confined  and  the  pore 
pressure  changes  might  reach  greater  depths  than  those 
predicted  by  these  investigations. 

5.  Porous  media  modelling  confirms  that  the  shallow 
activity  must  be  highly  related  to  the  pumping  of  the 
water . 

6.  The  shallow  activity  in  the  Valley  of  Mexico  seems  to 
have  no  relation  to  regional  stresses  that  could  derive 
from  intraplate  processes  since  the  location  of  unstable 
areas. was  not  modified  by  increments  in  the 

compress ional  regional  stress. 

The  results  obtained  here  for  the  Valley  of  Mexico  are 
in  excellent  agreement  with  the  hypothesis  that  diminishing 
the  amount  of  fluids  in  a  geologic  environment  can  increase 
the  instability  of  areas  prone  to  thrust  faulting,  and  may 


. 
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act  as  a  triggering  mechanism  for  earthquakes. 
Theoretically,  this  induced  activity  should  appear  almost 
instantaneously  after  the  extraction  of  fluids  takes  place, 
and  this  seems  to  be  the  case  in  this  area. 


4.  Thermo-elastic  Behaviour  of  the  Subducting  Lithosphere, 
and  its  Relation  to  Long  Range  Precursors. 

The  concept  of  instability  is  such  that  high  risk  areas 
will  appear  along  the  discontinuities  of  the  models  (Chapter 
1).  A  more  complex  problem  in  stability  studies  is  the 
distribution  of  instability  at  subducting  zones,  since  at 
the  Middle  America  trench  the  Cocos  plate  subducts  beneath 
the  other  two  plates  in  this  system.  Subduction  at  both 
sides  of  the  triple  junction  occurs  at  different  angles. 

In  this  chapter  different  subduction  models  are 
analysed  using  the  probability  function  described  in  section 
2.6.  Earthquakes  in  the  models  are  simulated  as  stress  drops 
at  the  focal  location.  The  analysis  of  several  different 
models  allows  the  suggestion  of  a  link  between  the 
characteristics  of  the  model  and  the  range  at  which  a 
seismic  event  can  affect  the  stability  of  other  parts  of  the 
model . 

Long  range  precursors  have  been  observed  along  some 
parts  of  the  seismic  belts  of  the  earth.  However,  the  range 
at  which  the  anomalous  seismicity  can  affect  future  activity 
varies  from  region  to  region.  The  state  of  stress  on  the 
lithosphere  also  has  strong  variations  through  the  different 
geographical  locations.  Here  the  possible  relationships  that 
exist  between  the  nature  of  the  material  that  constitutes 
the  lithosphere  and  the  occurrence  of  long  range  precursors 
will  be  explored. 
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The  probability  function  defined  in  section  2.6,  is 
applied  to  determine  the  risk  of  seismic  failure  in 
different  models  of  subduction  zones,  with  the  intention  of 
determining  what  consequences  the  occurrence  of  a  large 
earthquake  has  on  the  overall  probability  of  failure  in  the 
model . 

The  results  are  similar  to  the  seismicity  actually 
observed  in  this  kind  of  system.  They  strongly  suggest  that 
the  thermally  induced  stresses,  the  coupling  between  the  two 
converging  plates  and  the  angle  of  subsidence  may  play  a 
very  important  role  in  the  range  at  which  long  range 
precursors  are  possible.  The  range  of  these  seismic  patterns 
will  vary  widely  for  different  trenches. 

4 . 1  Introduction . 

The  determination  of  the  state  of  stress  in  the  earth's 
lithosphere  involves  a  number  of  plausible  assumptions;  One 
of  these  is  that  the  deep  continental  crust  is  less  dense 
than  deep  oceanic  crust.  Another  is  the  notion  that  the 
oceanic  lithosphere  increases  its  rigidity  and  its  thickness 
as  it  moves  away  from  the  ridge  at  which  it  was  generated 
(Leeds  et  al.,1974  and  Forsyth,  1979).  It  is  also  generally 
assumed  that  the  upper  part  of  the  plate  behaves  elastically 
and  only  when  its  temperature  increases  with  depth  does  it 
start  to  act  in  a  viscoelastic  way  (Solomon,  1980;  Wyss  1977 
and  others ) . 
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This  viscoelastic  behaviour  has  been  discussed  by 
Savage  and  Prescott  (1978),  who  pointed  out  that  the  short 
term  response  of  a  viscoelastic  material  will  satisfy  the 
requirements  for  the  upper  part  of  the  mantle.  The  short 
term  response  to  faulting  will  be  elastic.  He  also  concluded 
that  appreciable  stress  changes  may  be  produced  at  distances 
of  several  times  the  thickness  of  the  lithosphere  away  from 
the  source.  These  stress  alterations  are  relaxed  by  the 
material  flow,  so  that  the  response  of  the  viscoelastic 
material  is  an  instantaneous  pulse  followed  by  a  monotone 
relaxation . 

The  regularity  observed  in  earthquake  recurrence  in 
Japan  has  lead  Shimazaki  and  Nakata  (1980)  to  the  suggestion 
that  the  larger  an  earthquake  is  the  longer  the  following 
quiet  period.  This  allowed  them  to  propose  several 
recurrence  models  for  large  earthquakes.  To  do  this  they 
assumed  a  constant  accumulation  rate  of  tectonic  stress.  A 

similar  recurrence  pattern  has  been  observed  for  small 

* 

earthquakes  in  California  (Bufe,  1977).  The  constant  rate  of 
stress  accumulation  scheme  has  been  used  also  by  Newman  and 
Knopoff  (1982).  This  is  an  example  of  long  range  interaction 
in  time  in  an  earthquake  catalog.  Long  range  interaction  in 
space  may  also  exist. 

It  has  been  suggested  that  earthquakes  sometimes  affect 
the  seismicity  pattern  of  areas  at  great  distances  and  that 
some  large  earthquakes  seem  to  have  been  preceded  by  distant 
premonitory  seismicity  clusters.  Several  of  these  seismicity 
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patterns  for  long  range  precursors  have  been  proposed  by 
Keilis-Borok  et  al.  (1980  a,b,c).  These  patterns  consist  of 
an  abnormal  clustering  of  seismic  events  in  time,  energy  and 
space.  Pattern  Bursts  of  Aftershocks  consist  of  a  medium 
magnitude  earthquake  with  an  anomalous  sequence  of 
aftershocks.  Pattern  Swarms  consist  of  a  group  of  moderate 
events  concentrated  in  space  and  time.  Pattern  Sigma  roughly 
consist  of  an  increase  of  energy  released  in  a  sliding  time 
window  (Lamoreaux  et  al.,1983). 

These  patterns  have  been  applied  in  algorithm  form  with 
some  success  in  the  Central  Asia,  Eastern  Mediterranean, 
Pamir  and  Assam  regions  (Keilis-Borok  and  Malinovskaya, 

1964);  in  Italy  (Gasperini  et  al,  1978);  in  Anatolia  and 
surrounding  regions  (Keilis-Borok  and  Rotvain,  1979);  in  the 
Lake  Jocassee,  South  Carolina  region  (Sauber  and  Talwani, 
1980);  in  Southern  California  (Keilis-Borok  et  al.  1980  b); 
in  Tibet  and  the  Himalayas  (Keilis-Borok  et  al  1980  c);  and 
in  the  North  Amer ica-Cocos-Car ibbean  triple  junction  region 
(Lamoreaux  et  al,  1983). 

Premonitory  seismicity  patterns  are  regional  in  nature 
and  they  do  not  indicate  the  exact  location  or  the  time  of 
occurrence  of  a  major  earthquake.  The  results  must  be 
interpreted  as  an  increase  in  the  probability  within  a 
region  inside  a  time  window.  Physical  explanations  for  long 
range  precursors  are  controversial.  Barenblatt  et  al.  (1984) 
suggest  that  the  shear  stress  in  the  whole  area  is  large  and 
there  exist  local  areas  with  higher  concentrations.  They 
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claim  that  there  is  no  reason  to  expect  that  these 
concentrations  are  confined  to  the  vicinity  of  the 
earthquakes.  It  is  suggested  here  that  the  actual  physical 
properties  of  the  upper  part  of  the  lithosphere  may  very 
well  control  the  range  of  the  long  range  precursors. 

Newman  and  Knopoff  (1982),  on  the  other  hand,  pay 
little  attention  to  the  distribution  of  physical  properties 
and  suggest  that  the  presence  of  long  range  precursors  is 
linked  to  a  dynamical  process.  They  point  out  that  the 
fusion  of  small  cracks  into  larger  ones  can  produce  a 
cascade  of  events  which  may  culminate  in  a  massive 
earthquake,  due  to  the  nonlinear  character  of  this  process. 
In  a  later  paper  they  propose  that  by  assuming  that  the 
stress  accumulates  locally  at  a  uniform  rate  only  to  be 
released  by  a  large  earthquake,  the  rate  at  which  little 
cracks  are  produced  can  be  related  to  the  tectonic  related 
stress . 

The  acceptance  of  these  observed  seismicity  patterns 
runs  into  two  major  difficulties;  it  is  not  always  possible 
to  detect  them  and  the  range  at  which  they  might  have  some 
influence  on  the  local  seismicity  of  other  areas  varies 
greatly.  The  hypotheses  suggested  by  Barenblatt  et  al. 

(1984)  and  by  Newman  and  Knopoff  (1982)  deal  with  the  first 
problem.  In  this  chapter  a  mechanism  that  allows  a  different 
range  for  seismic  events  of  the  same  kind,  depending  on  the 
regional  material  properties  in  the  vicinity  of  their  origin 
will  be  proposed. 
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Most  large  earthquakes  along  island  arcs  are  thought  to 
be  a  result  of  fault  motion  across  the  zone  of  contact 
between  the  two  lithospheric  plates  that  form  the  subducting 
zone.  The  dimensions  of  these  slip  zones  vary  widely;  a 
rupture  zone  may  extend  as  much  as  1000  km  along  an  arc. 
Rupture  zones  of  extraordinary  length  tend  to  occur  along 
arcs  where  the  two  slabs  of  lithosphere  interact  across  a 
broad  interface.  Along  those  arcs  where  there  appear  to  be  a 
narrow  zone  of  contact  the  maximum  dimension  of  rupture  is 
limited  to  150  km  or  less  (Kelleher  et  al.,  1974). 

Kanamori  (1977)  suggested  that  it  is  almost  certain 
that  the  degree  of  mechanical  coupling  between  the  two 
lithospheric  slabs  varies  among  different  island  arcs.  An 
important  feature  brought  up  in  this  study  is  that  the 
rupture  lengths  of  great  earthquakes  of  nearly  identical 
magnitude  vary  greatly  for  different  regions.  Relationships 
between  the  type  of  interface  and  the  extent  of  the  rupture 
have  also  been  suggested  by  Isacks  et  al.  (1968),  and  for 
the  Ci rcum-Pac i f ic  belt  by  Kanamori  (1971). 

Kelleher  et  al.  (1974)  conclude,  that  the  maximum 
seismic  slip  that  can  occur  in  a  particular  trench  due  to  a 
large  earthquake  is  strongly  influenced  by  the  geometry  of 
the  interface  zone  and  in  particular  by  the  width  of  the 
plate  interfaces.  They  analyse  observations  from  several 
subducting  zones  and  from  their  study  the  following  remarks 
are  important  for  the  purposes  of  this  chapter. 

1.  Along  the  Middle-Amer ica  subduction  zone  where  the 
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extent  of  shallow  earthquakes  (<  70  km)  is  relatively 
narrow,  the  typical  observed  ruptures  are* of  70  to  150 
km  long. 

2.  Ruptures  of  extraordinary  length  have  occurred  along  arc 
segments  where  there  is  evidence  of  a  broad  interface 
between  both  slabs  of  lithosphere.  For  the  1957 
Aleutian,  the  1964  Alaskan  and  the  1960  Chilean 
earthquakes  ruptures  fault  length  extended  for  more  than 
800  km.  For  the  1952  Kamchatka  and  the  1965  Rat  Island 
earthquakes  ruptures  length  extended  to  more  than  400 

km . 

3.  Near  each  of  these  large  rupture  zones  along  island 
arcs,  the  line  of  volcanoes  approaches  the  surface 
projection  of  the  contour  line  for  shallow  earthquakes 
(<  70  km).  By  contrast,  these  two  lines  are  clearly 
separated  along  arc  segments  characterized  by  moderate 
rupture  zones. 

The  consistency  of  these  observations  suggest  that 
there  should  be  a  physical  relationship  between  the  geometry 
of  the  subducting  zone,  the  location  of  active  volcanoes  and 
the  extent  of  the  ruptures  due  to  large  shallow  earthquakes. 

On  the  basis  of  multiple  sei sinological  observations 
Kanamori  (1977),  proposed  that  a  parameter  that  can  be 
responsible  for  the  difference  in  rupture  lengths  might  very 
well  be  the  degree  of  coupling  between  the  two  lithospheric 
slabs  involved  in  the  subduction.  He  suggested  a  plate 
coupling  and  decoupling  model  as  an  attempt  to  understand 
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the  physical  mechanisms  in  the  various  types  of  subducting 
zones  (figure  21).  This  model  is  an  extension  of  the  one 
proposed  by  Kanamori  (1971). 

Figure  21  (Kanamori,  1977;  figure  8)  shows  the 
different  kinds  of  possible  subducting  mechanisms  that 
result  from  coupling  considerations.  Figure  21a  represents  a 
strong  coupling  between  the  oceanic  and  continental 
lithospheres  which  results  in  great  earthquakes  and 
break-off  of  the  subducting  lithosphere  at  shallow  depths. 

In  21b  partial  decoupling  leads  to  smaller  earthquakes  and 
continuous  subduction  while  in  21c  further  decoupling  might 
lead  to  aseismic  events  or  intraplate  tensional  events.  In 
21d  the  sinking  of  the  plate  increases  the  angle  of 
subduction  and  results  in  retreating  subduction  and  the 
formation  of  a  new  thin  layer  of  lithosphere. 

4.2  The  Method. 

The  initial  state  of  stress  is  of  great  importance  in 
the  determination  of  the  risk  of  reaching  failure,  therefore 
the  shear  stress  concentrations  that  may  be  present  due  to 
faults  or  other  structures  are  areas  with  high  potential  for 
failure.  The  stress  build-up  rate  as  discussed  by  Shimazaki 
and  Nakata  (1980)  and  by  Newman  and  Knopoff  (1982)  justifies 
the  process  of  modelling  a  region  before  and  after  the 
occurrence  of  a  stress  drop  due  to  an  earthquake  as  a  static 
problem. 
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Figure  21....  Different  types  of  subducting  zones  resulting 
from  plate  coupling  considerations  (modified  from  Kanamori, 
1977). 


Such  modelling  requires  computation,  almost  inevitably 
by  finite  element  methods.  The  finite  element  technique  is 
one  of  the  most  closely  scrutinized  methods  in  stress 
calculation,  since,  as  in  all  quantitative  analyses,  the 
theory  on  which  it  is  based  makes  approximations  and  several 
different  approaches  can  be  taken  for  a  particular  problem. 
Also  the  application  of  realistic  boundary  conditions  is 
often  a  dificult  task.  Only  through  experience  can  the  best 
choice  be  made.  As  a  result,  conclusions  derived  from 
numerical  techniques  are  often  viewed  with  reserve.  It  is 
obviously  desirable  to  solve  the  same  problem  by  other 
methods  to  serve  as  an  independent  check. 

The  temperature  distribution  used  for  both  of  the 
plates  involved  in  the  subducting  mechanism  here  is  in 
accordance  with  that  suggested  by  McKenzie  (1969).  Where  the 
temperatures  inside  the  subducting  slab  are  calculated  from 

3  2  T '  3T '  3  2  T ? 

- - -2R- - + - =0  (4.1) 
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where  T=T1T',x=lx'  and  z=lz'.  The  solution  of  this  equation 
is  of  the  form 

T’ =  1 +2Zn  ( -  1  )  n  •  exp{R- (R2 +n  27r 2  )  1  7  2  }x  '  sin  n7rz'/n7r  (4.2) 

In  this  formulation  R,  the  Reynolds  number,  does  not  have  y 
dependence  because  the  slab  is  assumed  a  two  dimensional 
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structure . 

The  values  taken  for  the  parameters  in  this  chapter 

are : 

Specific  heat  Co  =  0.25  cal  g",cC'' 

Thermal  expansion  coefficient  a=4x 1 0 " 5 cC' ' , 

Temperature  T • =800 °C 

Thermal  conductivity  k=.001  cal  cm' 1 °C' 1  sec ' ‘ 

Width  of  the  slab  1=40  km. 


The  resulting  t 
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the  nodes  of  the  model. 

It  is  important  to  start  the  analysis 
model  in  order  to  obtain  a  general  idea  of 
changes  as  the  model  becomes  more  ccmplicat 
modelling  has  been  proven  useful  in  several 
problems  (i.e.  Solomon,  1980;  Hanks  and  Ral 
in  particular  in  seismic  stability  studies 
and  Nv land ,  1  983  a ,  b ) .  Nevertheless  elastic 
be  representative  of  the  upper  few  tens  of 
earth's  crust,  where  the  temperature  of  the 
material  is  beneath  300 c C .  Below  this  tempe 
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thermally  induced  stress  is  small  and  the  tectonics  is 
describable  by  elastic  stresses  (Turcotte,  1974). 

The  lithosphere  can  be  thought  of  as  an  elastic  layer 
over  a  high  viscosity  layer.  Beneath  the  lithosphere  there 
exists  a  'low  strength'  zone  known  as  the  asthenosphere , 
which  by  definition  is  not  capable  of  resisting  great  loads 
for  a  very  long  period  of  time.  This  means  that  at  a  given 
depth,  above  the  asthenosphere,  the  stress  due  to  permanent 
vertical  loads  cannot  produce  vertical  displacements  after  a 
certain  relaxation  time.  It  is  unlikely  that  the  boundary 
between  the  lithosphere  and  the  asthenosphere  is  sharp;  the 
lower  part  of  the  lithosphere  should  have  similar  physical 
properties  to  those  of  the  upper  part  of  the  asthenosphere. 
The  material  in  the  upper  layer  of  the  plate  should  become 
gradually  less  viscous  as  depth  decreases. 

Since  oceanic  lithosphere  cools  as  it  moves  away  from 
the  ridges  (Forsyth,  1979)  its  strength  and  density 
increases.  Strength  and  density  are  the  main  difference 
between  continental  and  oceanic  plates  which  is  primarly  due 
to  the  physical  characteristics  of  their  materials.  Strength 
is  the  maximum  differential  stress  (amx-amin)  that  can  be 
applied  to  the  material  without  reaching  failure;  it 
increases  with  increasing  confining  pressure.  It  has  been 
assumed  here  that  the  oceanic  lithosphere  has  basaltic 
physical  properties,  while  continental  plates  have 
rhyolitic . 
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The  most  important  force  acting  in  a  boundary  zone 

between  two  plates  with  converging  relative  velocities  is 

the  gravitational  pull  of  the  downgoing  slab  (Chappie  and 

Tullis,  1977).  This  force  has  been  proposed  as  the  probable 

\ 

driving  mechanism  of  plate  tectonics  (e.g.  Press,  1972; 
McKenzie,  1972;  Richter  1977).  The  frictional  forces  acting 
on  plate  boundaries  are  very  small  compared  with  the 
gravitational  pull,  the  viscous  drag  of  the  mantle,  and  the 
ridge  push  force  (Chappie  and  Tullis,  1977).  However,  it  is 
of  great  importance  to  include  these  frictional  forces  in 
stability  studies,  since  it  has  been  proven  in  our  previous 
studies  (Ur ibe-Carva jal  and  Nyland,  1982,1983)  that  small 
shear  stress  changes  can  affect  the  overall  stability 
distribution.  The  shear  acting  at  the  boundaries  of  all  our 
models  is  in  accordance  with  the  one  determined  by  McKenzie 
(1969). 

The  other  force  that  must  be  included  in  any  subduction 
zone  model  is  the  viscous  drag  of  the  mantle;  it  is  the  only 
one  with  the  capability  of  balancing  the  action  of  the 
slab-pull  (Richardson,  1978),  so  that  plate  velocities  are 
restricted  to  few  centimeters  per  year.  Such  velocities  are 
indicated  by  Chase  (1972),  Minster  et  al.  (1974),  Sacks 
(1983)  and  others,  for  all  the  plates. 

The  resistive  forces  that  oppose  the  movement  of  the 
plates  are  due  to  differential  movements  at  the  edges  of  the 
plates  and  viscous  counter  flow  of  the  asthenosphere .  The 
resistive  forces  at  the  interplate  boundaries  are  difficult 
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to  estimate.  However,  earthquakes  clearly  show  that  these 
forces  are  not  sufficient  to  stop  the  movement  of  the 
plates.  This  has  been  recognized  by  Jacoby  (1970), 

Richardson  (1978)  and  many  others. 

The  assumption  that  there  is  an  elastic  lithosphere 
drifting  over  a  viscous  asthenosphere  implies  the  existence 
of  a  drag  force  opposing  this  movement.  In  order  to  model 
this  force  it  is  necessary  to  assume  a  velocity  distribution 
in  the  asthenosphere  that  depends  on  its  flow.  If  the 
elastic  plate  is  drifting  over  a  passive  viscous  layer  it  is 
dragged  over  a  base  which  is  at  rest. 

The  drag  stress  per  unit  length  (Richardson,  1978)  is 

T  =  riV/h  (4.3) 

where  V  is  the  velocity  of  the  plate,  77  and  h  the  viscocity 
and  the  thickness  of  the  asthenosphere.  The  term 
asthenosphere  applies  to  the  viscous  layer  that  underlies 
the  elastic  plate,  at  the  bottom  of  this  plate  no 
displacement  is  possible,  and  the  drag  force  is  then  a 
function  of  its  thickness  ih) .  Following  Richardson  (1978), 
we  take  a  velocity  V=10cm/year,  h=^  3x107cm  and  77=  1  0  2  0  poise, 
to  obtain  r=1bar. 

The  integration  of  this  force  over  the  bottom  surface 
of  a  plate  gives  values  that  can  equlibrate  the  stresses  due 
to  the  slab  pull  and  the  ridge  push.  The  frictional  forces 
due  to  the  plate  boundaries  and  transform  faults  are  much 
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smaller  and  represent  a  secondary  contribution  to  the 
overall  quasi-equilibrium  of  the  plates.  All  these  forces 
have  been  calculated  for  all  the  major  plates  by  Chappie  and 
Tullis  (1977).  They  demonstrate  that  only  a  few  forces  are 
responsible  for  the  observed  velocities  of  the  plates.  These 
are  the  ones  that  have  been  taken  into  account  in  our  model 
calculations . 

In  this  chapter  a  geodynamical  study  of  a  subduction 
zone  is  presented.  This  is  done  using  the  ADINA  finite 
element  program  (Bathe,  1978).  The  use  of  models  is  only 
justified  to  the  extent  that  the  model  is  a  reasonable 
facsimile  to  reality,  and  in  order  to  be  able  to  model 
stresses  in  a  subduction  zone  several  assumptions  have  to  be 
made.  The  stress-strain  relationships  that  accurately 
represent  the  stresses  at  depth  are  uncertain.  Nevertheless, 
for  the  use  of  stability  functions  it  is  sufficient  to  have 
an  approximate  model  since  they  are  affected  only  locally  by 
small  variations  to  the  model  (Ur ibe-Carva jal  and  Nyland, 
1982, 1983) . 

In  the  present  case  a  continuum  in  which  material 
property  changes  can  be  modelled  by  finite  elements  is  dealt 
with,  and  the  boundary  conditions  can  be  specified  at  the 
required  locations.  Therefore,  the  internal  properties  of 
the  model  can  be  varied  and  the  variations  of  stability 
(derived  from  the  stresses  determined  by  ADINA)  that  result 
can  be  interpreted.  Finite  element  models  of  subducting 
zones  have  been  studied  previously  by  others.  Bird  (1978), 
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applies  the  technique  to  calculate  the  shear  stress 
distribution  within  a  subduction  system  involving  a  slab 
dipping  at  45°.  The  greatest  difference  between  Bird's  model 
and  the  one  shown  here  is  the  boundary  condition  at  the 
upper  part  of  the  model;  he  considers  the  plates  to  be  fixed 
at  the  surface,  that  is,  the  model  does  not  allow  for 
gravitational  deformation  of  the  plates. 

The  present  model  assumes  that  the  system  is  not 
hanging  from  the  surface,  but  supported  by  a  surface  at  its 
bottom.  The  same  kind  of  boundary  condition  was  used  in  the 
analysis  of  the  Valley  of  Mexico  in  chapter  3  of  this 
Thesis . 

The  model  presented  here  consists  of  an  oceanic  layer 
of  40  km  of  thickness  underthrusting  a  80  km  thick 
continental  layer.  In  the  lower  part  of  these  plates  the 
temperatures  are  well  above  the  300°C  value  (i.e.  McKenzie, 
1969?  Griggs,  1972;  Minear  and  Toksoz  1970),  accepted  as  the 

limit  for  elastic  behaviour.  The  thermal  dependence  of  the 

✓ 

stresses  plays  an  increasingly  important  role  in  the  overall 
state  of  stress  as  depth  becomes  greater.  This  can  be 
modelled  using  the  ADINA  option  for  thermo-elastic 
modelling.  In  essence,  the  variation  with  temperature  of  the 
elastic  parameters  is  specified. 

The  thermal  dependence  of  the  Young's  modulus  (E)  is 
taken  to  be 

E=Em  j  n  +E0e " x 


(4.4) 
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where  E0=E  at  the  surface  of  the  earth  minus  Em s n ,  x=T2/C,  T 
is  the  temperature  and  C=600.  This  equation  describes  the 
exponential  non-linear  temperature  versus  depth  behaviour  of 
the  upper  mantle  suggested  by  Kirby  (1977)  and  others.  It 
should  be  noted  that  as  depth  increases  the  temperature  of 
the  material  is  also  increasing.  At  small  temperatures  E 
remains  virtually  unchanged  up  to  a  given  temperature  T,, 
where  the  material  starts  to  increase  its  viscous  behaviour 
rapidly  until  a  temperature  T2  is  reached.  Then  the  value  of 
E  becomes  asymptotically  equal  to  a  lower  limit  value  Em  ■,  n . 
It  is  our  claim  that  this  Em j n  is  that  one  representative  of 
the  value  of  viscosity  of  the  upper  asthenosphere  stated 
earlier  as  77. 

The  exponential  decay  of  E  with  temperature  has  been 
suggested  previously  in  geophysics.  Vening  Meinesz  (1964) 
pointed  out  that  at  a  certain  temperature  the  elastic  limit 
indeed  diminishes  and  that  for  some  substances  like  Mg2Si04, 
which  is  assumed  to  be  an  important  constituent  of  the 
mantle  at  a  given  temperature,  this  elastic  limit  does  not 
continue  to  decay. 

Theoretically,  Poisson  ratio  v  has  an  upper  limit  of 
0.5,  and  for  most  crustal  rocks  v  is  around  0.25  and  it  is  a 
common  assumption  in  geophysical  modelling  to  use  this  0.25 
value.  Therefore,  I  have  adopted  this  value  for  surface 
rocks  and  have  interpolated  to  p=0.3  for  the  rocks  at  the 
bottom  of  the  plate.  The  value  of  the  thermal  conductivity, 
a,  has  been  taken  to  be  constant  (McKenzie,  1970)  throughout 
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the  whole  exercise.  The  values  of  the  parameters  have  been 
varied  (table  1)  to  simulate  the  dependence  of  temperature 
variations  on  the  materials  modelled. 

There  is  an  apparent  contradiction  between  the  values 
in  table  1  and  the  general  notion  that  seismic  velocities 
increase  as  depth  increases  (Stacey,  1977,  figure  6.18),  and 
the  increasing  Young's  modulus  E  with  depth  used  by  Turcotte 
(1974).  This  work  focusses  on  the  behaviour  of  the  lower 
part  of  the  plate  from  a  quasi-static  or  isothermal  point  of 
view.  We  are  not  dealing  with  the  dynamical  or  adiabatic 
properties  of  the  material.  Then  as  depth  increases  E  should 
become  smaller  because  the  material  is  becoming  more  ductile 
as  temperature  increases,  even  if  the  moduli  are  not 
explicitly  temperature  dependent. 

The  difference  between  the  adiabatic  (dynamic)  and 
isothermal  (static)  parameters  comes  from  allowing 
temperature  changes  within  the  model  of  the  earth.  The 
relation  between  them- can  be  found  directly  from 

Ed=E/( 1 -ETa  2 /9C0 ) 

v  d  =  ( j'+ETa  2/9C o  )  / (  1  ~ETa  2 /9C 0  )  (4.8) 

where  C0  is  the  specific  heat  per  unit  volume  at  constant 
pressure.  Ed  and  v  d  are  the  dynamic  Young's  modulus  and 
Poisson's  ratio  (Landau  and  Lifshitz,  1970).  Nevertheless, 
this  kind  of  differentiation  between  dynamic  and  static 
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DENSITY 

3. 

TEMP. 

20 

100 

YOUNG 

.  7E  12 

. 68E  12 

POISSON 

.25 

.25 

ALPHA 

4.E-5 

4.E-5 

DENSITY 

2.7 

TEMP. 

20 

100 

YOUNG 

. 6E+ 1 2 

. 58E 1 2 

POISSON 

.25 

.25 

ALPHA 

4.E-5 

4.E-5 

SUBDUCTING  PLATE 


200 

400 

600 

63E12  . 

46E12 

.  IE  12 

.26 

.265 

.27 

4.E-5 

4.E-5 

4.E-5 

UPPER 

PLATE 

200 

400 

600 

54E12  . 

40E12 

. 24E  1 2 

.26 

.265 

.27 

4.E-5 

4.E-5 

4.E-5 

800 

1000 

1250 

.  7E  1 1 

. 57E 1 1 

.  3E  1 1 

.28 

.29 

.3 

4.E-5 

4.E-5 

4.E-5 

800 

1000 

1250 

12E12 

. 49E 1 1 

.  2E  1 1 

.28 

.29 

.3 

4.E-5 

4.E-5 

4.E-5 

Table  1....  Shows  the  dependence  of  the  material  properties 
on  temperature  of  the  thermo-elastic  models  studied. 

Density  is  given  in  g/cm3 ,  Young  modulus  in  dyne/cm2, 
Temperature  in  °C,  and  a  in  0C-1 
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parameters  only  makes  their  values  vary  by  a  small 
percentage.  An  additional  statement  is  made  however,  by 
attempting  to  include  softening  due  to  increasing 
temperature . 

All  the  forces  involved  in  the  models  are  placed  at  the 
node  nearest  to  the  actual  point  of  application,  as  is  done 
in  other  geophysical  finite  element  modelling  exercises 
(Richardson,  1978;  Richardson,  Solomon  and  Sleep,  1979);  The 
occurrence  of  an  earthquake  is  modelled  by  introducing  a 
sudden  displacement  or  stress  drop  in  the  model  at  the 
location  of  the  hypocenter.  This  kind  of  earthquake 
simulation  is  usual  in  finite  element  analysis  (Richardson, 
1974).  The  stress  drop  due  to  an  intermediate  depth 
intraplate  earthquake  with  focus  in  a  subduction  zone  is  of 
the  order  of  100  to  200  bars  (Chinnery,  1964).  Recently  this 
value  is  now  believed  to  be  arround  40  bars. 

4.3  The  subduction  zone  models. 

Several  types  of  island  arcs  have  been  discussed  and 
and  summarized  in  figure  21.  Geometrical  changes  in  the 
model  and  the  variation  of  the  forces  involved  in  the  model 
should  generate  changes  in  the  overall  distribution  of 
seismic  instability.  If  the  seismic  behaviour  of  the 
different  types  of  subduction  zones  depends  mostly  on  these 
parameters,  it  is  then  necessary  to  study  how  these  changes 
are  reflected  in  the  probability  distribution. 
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This  kind  of  geometrical  dependence  of  any  risk 
*  function  must  result  from  the  analysis  of  several  models. 

The  models  considered  here  will  include 

1.  Material  property  changes,  from  elastic  to 
thermo-elastic  stress-strain  relationships. 

2.  Changes  in  the  angle  of  subduction  and  the  value  of  the 
frictional  forces.  This  will  account  for  the  coupling 
and  partial  uncoupling  effects  (Kanamori,  1977). 

3.  The  length  of  the  subducting  slab  has  also  been  changed 
in  order  to  model  slabs  of  different  ages,  which  also 
involves  changes  of  the  gravitational  pull  exerted  by 
the  downgoing  slab. 

The  element  distribution  (figure  22)  used  for  the  model 
studies  shown  in  figures  23,  24  and  25  is  the  same.  However, 
the  ratio  between  the  horizontal  and  vertical  total  lengths 
varies  to  model  slabs  subducting  at  different  angles. 

The  initial  model  is  one  plate  underthrusting  a  second 
plate  at  an  angle  of  45°  (figure  23).  A  linear  elastic  model 
of  this  kind  of  subduction  zone  and  the  boundary  conditions 
and  restrictions  mentioned  produce  the  probability 
distribution  given  in  figure  23a.  The  higher  risk  of  failure 
seems  to  be  concentrated  in  specific  locations,  at  which  the 
analysis  will  be  focused. 

The  high  risk  area  is  located  in  the  intraplate 
boundary  from  a  depth  of  about  50  km.  This  area  has  a 
probability  of  failure  of  more  than  0.75  with  respect  to  the 
other  parts  of  the  model.  The  highest  probability  of  failure 
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Figure  22 ... . 


Nodal  distribution  used  in  the  initial  models 
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Figure  23....  Elastic  modelling  of  a  subducting  zone,  a,  b 
are  the  probability  distributions  before  and  after  the 
occurrence  of  an  earthquake  at  *. 
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occurs  at  the  location  *  on  this  figure.  It  seems  from  this 
simple  model  that  there  tends  to  be  a  band  of  instability  at 
depth  along  the  boundary,  which  strongly  suggests  that  the 
band  should  be  present  for  models  containing  longer  (older) 
subducting  slabs. 

Figures  23a  and  23b  represent  the  same  model  and  the 
only  difference  is  that  in  figure  23b  a  seismic  event  has 
been  introduced  in  the  location  *  which  is  inside  the  high 
risk  valued  zone  of  figure  23a.  The  seismic  event  has  been 
simulated  by  introducing  a  stress  drop  of  100  bars  at  the 
desired  point  of  the  model.  This  value  corresponds  to  that 
calculated  for  intraplate  earthquakes  of  magnitude  7  ,  with 
^=3x10’  2  dynes/cm2,  77=. 75,  D=5  m,  and  w=100  km,  using  the 
Kanamori  and  Anderson  (1975)  formula  for  a  strike-slip 
mechanism: 

Aa  =  2juD/w7r  (4.9) 

where  Ao  is  the  stress  drop  due  to  the  earthquake,  ju  is  the 
rigidity  of  the  material,  D  is  the  average  offset  and  w  is 
the  rupture  width. 

The  probability  distribution  throughout  the  model  has 
been  changed.  After  the  simulation  of  the  earthquake  the 
area  of  the  hypocenter  has  become  more  stable  since  stress 
has  been  released  there.  The  high  valued  zone  seems  to  have 
been  broken  apart  and  expanded  in  a  radial  direction.  The 
two  lower  high  valued  zones  might  correspond  to  a  double 
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Benioff  zone  similar  to  that  observed  in  Japan  (Uyeda, 

1977),  and  in  the  Central  Aleutian  Arc  (Engdahl  and  Scholz, 
1977).  The  minimum  stability  zone  has  moved  nearer  to  the 
inside  of  the  slab.  This  result  implies  that  the  occurrence 
of  an  earthquake  in  this  kind  of  model  indicates,  by 
changing  the  whole  probability  distribution,  where  the  next 
failure  is  likely  to  take  place.  It  is  a  comparative  scheme 
and  shows  that  the  introduction  of  stress  changes  may  not 
directly  affect  particular  stress  components  in  the  distant 
zones  significantly,  but  it  can  produce  changes  in  the 
probability  of  failure  with  respect  to  the  model  as  a  whole. 

The  overall  stability  of  a  region  is  changed  by  the 
occurrence  of  a  large  earthquake;  stress  accumulates 
constantly  and  earthquakes  are  required  to  release  this 
stress.  Therefore,  changes  in  the  probability  function  may 
very  well  represent  new  states  of  quasi-equilibrium  for  the 
whole  system,  since  it  is  defined  as  the  probability  of 
failure  of  a  given  location  with  respect  to  the  whole  model. 

A  smaller  stress  drop  was  tried  to  see  how  far  the 
changes  due  to  such  a  seismic  event  can  be  seen.  A  stress 
drop  of  30  bars  was  then  modelled  and  the  general  results 
were  the  same  -  a  general  reduction  in  probability 
throughout  the  model  and  a  small  migration  of  the  high 
valued  area  towards  the  inside  of  the  slab. 

From  the  discussion  in  the  introductory  part  it  can  be 
seen  that  the  elastic  model  of  the  whole  plate  is 
unrealistic.  An  improvement  is  the  analysis  of  a 
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thermo-elastic  model  (figure  24)  where  the  same  approach  as 
in  the  elastic  case  was  taken.  Figure  24a  shows  the 
probability  function  distribution  for  a  simple 
thermo-elastic  model  and  figure  24b  shows  the  same  model 
including  a  stress  drop  of  100  bars  at  the  location  shown. 

Figure  24a  shows  that  in  the  plane  part  of  the 
underthrust ing  plate  there  is  a  high  probability  of  failure. 
This  pattern  remains  for  the  upper  part  of  the  plate  until 
depths  at  which  the  material  has  almost  reached  its  plastic 
state.  High  valued  zones  also  appear  in  the  same  region  as 
that  in  the  elastic  model  and  in  a  narrow  zone  which  runs 
inside  the  subducting  slab  approximately  parallel  to  its 
dip.  In  this  figure  again  the  double  Benioff  zone  seems  to 
be  present,  and  the  area  with  major  probability  is  located 
in  the  intraplate  boundary  at  approximately  60  km  deep.  As 
in  the  elastic  model  after  an  explosion  of  100  bars,  the 
general  features  of  the  probabilistic  distribution  have  been 
altered;  the  probability  in  the  elastic  part  of  the  slab  has 
suffered  the  greatest  change  and  the  stability  in  the 
viscous  parts  remains  greater  than  in  the  elastic  layers  of 
the  plates.  The  viscous  parts  have  not  been  affected  at  all 
except  in  the  immediate  neighbourhood  where  the  explosion 
took  place.  Another  similarity  observed  between  this  and  the 
previous  model  is  the  presence  of  some  kind  of  ring  of  high 
risk  areas  surrounding  the  stress  drop  and  that  aftershocks 
on  both  plates  are  likely  to  occur. 
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Figure  24....  Thermo-elastic  model 
45°.  a  and  b  the  same  as  in  figure 


of  a  subducting  zone  at 
23. 
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From  this  comparison  it  can  be  seen  that  the  upper  part 
of  both  plates  becomes  relatively  more  unstable  in  the 
thermal  analysis.  This  indicates  that  the  viscous  part  has 
actually  become  more  stable  than  in  the  elastic  case,  since 
the  parameters  of  the  elastic  parts  remain  the  same.  After 
the  occurrence  of  the  stress  drop  the  "instability  ring” 
seems  to  have  propagated  further  away  from  the  focus  in  the 
elastic  case.  This  implies  that  in  the  viscous  material 
stress  propagation  will  have  a  larger  attenuation. 

This  result  is  for  a  plate  model  consisting  of  an 
elastic  layer  overlying  a  material  which  is  gradually 
becoming  viscous  with  a  thermal  dependence  for  Young's 
modulus  given  by  equation  (4.4).  It  strongly  suggests  that 
the  effect  that  a  given  stress  change  has  on  the  probability 
distribution  is  greatly  dependent  on  the  material  properties 
assumed.  It  presents  the  possibility  that  changes  in  this 
function  can  be  observed  at  distant  locations  from  the 
source  if  the  stress  change  can  be  transmitted  to  them 
through  elastic  materials,  i.e.  viscous  materials  should  act 
as  highly  attenuate  the  shear  stress  and  shear  energy 
transmission  within  the  kind  of  static  that  are  dealt  with 
here . 


The  next  step  is 
probability  function 
angle  of  subsidence, 
(chapter  3 ) ,  this  is 
how  geometry  affects 


to  determe  the  effects  on  the 
distribution  of  the  variation  of  the 
As  in  the  case  of  the  Valley  of  Mexico 
necessary  to  create  a  general  idea  of 
instability,  and  to  investigate  why  the 
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different  kinds  of  trench-arc  systems  described  in  figure  21 
have  different  observed  seismic  behaviour. 

Probability  distributions  were  calculated  for  two 
models  using  the  same  kinds  of  boundary  conditions  used  for 
figure  24,  the  temperature  dependence  of  the  material 
properties  and  the  temperature  distribution  of  the  models 
being  maintained  the  same. 

The  angle  at  which  the  slab  subducts  was  changed  to  10° 
(figure  25),  and  the  intraplate  frictional  forces  were  made 
3/2  of  that  used  in  the  previous  model.  This  is  an  effort  to 
model  the  high  coupling  which  is  expected  in  low  angle 
subduction  (Kanamori,  1977).  Figures  25a  and  25b  represent 
the  probability  distribution  before  and  after  the  modelling 
of  a  100  bar  stress  drop  at  location  *.  The  high  instability 
zone  (valued  at  0.6  for  illustration  purposes)  shows  a  more 
continuous  distribution  than  that  of  figure  24a,  and  the 
highest  valued  area  (0.8)  seems  to  be  locally  concentrated. 
Again,  most  of  the  unstable  area  coincides  with  the  elastic 
material  in  the  upper  part  of  the  lithosphere  and  the 
elastic  part  of  the  subducting  slab.  However,  when  the 
earthquake  has  been  introduced,  the  instability  seems  to 
have  moved  to  the  viscous  parts.  The  exception  to  this 
general  result  is  the  zone  above  the  hypocenter  which  may 
correspond  to  the  upper  part  of  the  instability  ring 
observed  in  figure  24b  . 

The  next  model  (figure  26)  is  the  result  of  changing 
the  angle  of  subduction  of  the  downgoing  slab  to  60°  and  the 
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Figure  25....  Probability 
slab  subducting  at  10°.  a 


distributions  for 
and  b  as  before. 
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trench  with  a 
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o?9Kre  Pr0babiiity  distributions  for  a  trench  with  a 

slab  subducting  at  60  .  a  and  b  as  before. 
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frictional  intraplate  force  to  1/2  of  that  of  figure  25. 

This  last  modification  is  made  to  simulate  the  partial 
uncoupling  between  the  two  plates,  described  earlier  for 
this  kind  of  model.  Figure  26a  (before  the  earthquake)  shows 
similar  characteristics  to  figure  24a,  the  elastic  layer  of 
the  lithosphere  and  a  region  in  the  intraplate  boundary 
having  the  highest  probability  of  reaching  failure.  The 
simulation  of  the  earthquake  (figure  26b)  shows  that 
although  the  area  surrounding  the  focus  has  become  the  most 
stable,  this  effect  seems  to  be  highly  localized.  The  upper 
part  of  both  plates  remains  approximately  at  the  same 
probability  as  in  figure  26a  .  The  lower  high  seems  to  have 
been  pushed  upwards. 

Older  slabs  can  be  modelled  by  applying  the  same  kind 
of  reasoning  to  larger  models.  The  results  of  using  a  model 
with  a  slab  dipping  at  45°  (figure  27),  and  one  with  the 
slab  at  60°  (figure  28)  are  presented.  Again  a  and  b 
represent  the  probability  distribution  before  and  after  the 
occurrence  of  a  large  event  at  location  *  which  is  the 
location  with  maximum  probability.  The  general  results  are 
similar  to  those  of  figures  25  and  26,  but  the  presence  of 
this  kind  of  high  probability  at  these  depths  requires  more 
discussion  since  at  these  depths  the  earthquake  mechanisms 
are  not  easily  visualized  as  sudden  failures  due  to  applied 
stresses  but  rather  to  volume  increases  due  to  phase  changes 
suffered  by  the  subducting  material. 
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Figure  27....  Probability  distributions  for  a  trench  with  a 
long  slab  subducting  at  45°.  a  and  b  as  before. 
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Figure  28....  Probability  distributions  for  a  trench  with  a 
long  slab  subducting  at  60°.  a  and  b  as  before. 
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4.4  Discussion  of  results. 

A  change  in  the  overall  distribution  of  the  probability 
function  values  does  not  mean  that  the  stresses  in  a 
particular  location  have  changed  substantially,  but  rather 
that  its  probability  of  reaching  failure  has  been  modified. 
The  degree  to  which  a  point  is  affected  depends  on  the 
physical  properties  of  the  material  in  which  it  is  located 
and  the  state  of  stress  in  which  it  was  previously. 

If  in  fact  the  shear  energy  is  trapped  in  the  elastic 
part  of  the  plate  then  there  exists  the  possibility  that  the 
thickness  of  this  elastic  layer  (H)  will  determine  the  range 
at_  which  the  effects  of  a  stress  drop  (like  the  earthquake) 
should  affect  the  probability  of  failure.  Then  it  is 
possible  that  the  reason  for  the  existence  of  the  so-called 
long  range  precursors  is  linked  directly  to  H.  It  is 
possible  that  the  stress  drop  due  to  an  earthquake  produces 
changes  at  greater  distances  than  those  predicted  by  simple 
theory . 

To  test  this  hypothesis  it  is  necessary  to  calculate 
the  solution  to  the  Boussinesq’s  problem  for  the  case  of  an 
elastic  layer  over  a  viscoelastic  half-space.  Solutions  to 
this  problem  have  been  known  since  the  begining  of  the 
century  (Michell,  1900). 

The  shear  energy  is  (Landau  and  Lifshitz,  1970,  pi  1  ) 


F=n ( e j  j -  1/3  6  j  j  e  k  k  )  2 


(4.10) 


n 
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where  e  ,  }  are  the  components  of  the  strain  tensor.  This 
exercise  gives  the  result  that  for  all  H  (within  reasonable 
limits)  the  total  shear  energy  in  the  elastic  plate  is  more 
than  0.9  of  the  total  shear  energy  of  the  whole  model.  This 
is  expected  from  our  finite  element  models. 

The  result  that  most  of  the  shear  energy  tends  to 
remain  within  the  elastic  part  of  the  plate  can  be  useful  to 
test  our  stability  results  obtained  by  finite  element 
calculations.  Both  problems  are  similar  and  can  be  solved  by 
making  the  same  assumptions,  using  the  same  boundary 
conditions  and  same  kind  of  material  elements.  A  finite 
element  calculation  of  the  same  Boussinesq's  problem  has 
been  prepared  by  assuming  a  two  layered  material.  The 
elastic  constants  chosen  are  E,=7x1012,  E2=7x1010,  v  ,  =  .25 
and  i>2=*3.  These  correspond  to  temperatures  of  20°C  and 
800°C  on  our  thermo-elastic  subduct  ion  model.  The  upper 
layer  has  a  thickness  of  40  km  and  the  bottom  one  90  km.  The 

shear  energy  in  the  upper  layer  is  again  more  than  90%  of 

* 

the  total  shear  energy  of  the  model.  This  result  adds 
credibility  to  the  previous  stability  calculations. 

The  concentration  of  shear  energy  in  the  elastic  layer 
is  a  common  feature  of  both  of  the  solutions  to  the 
Boussinesq's  problem.  It  is  therefore  possible  without 
losing  generality  to  approach  this  problem  using  the 
continuum  mechanics  thin  plate  theory  as  described  by  Love 
(1924),  Landau  and  Lifshitz  (1970)  and  many  other  textbooks. 


. 
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The  specific  problem  to  solve  is:  For  a  thin  plate 
determine  the  effects  of  an  applied  force  at  a  given 
distance  r;  then,  by  changing  H,  the  distance  r,  at  which 
the  same  effect  calculated  previously  can  be  determined  in 
such  a  way  that  the  variation  between  H  and  the  range  at 
which  a  fixed  observation  occurs  can  be  stated. 

The  deflection  of  a  circular  plate  with  clamped  edges 
due  to  the  application  of  a  force  f  in  its  center  is  given 
by 

$  =  3f  {  l-v2  )/27rEh3  [1/2  ( R2 -r  2  )  In  ( R/r  )  ]  (4.11) 

where  h=H/2  and  R  is  the  radius  of  the  plate. 

It  can  easily  be  demonstrated  that  if  $  is  to  remain 
constant  as  H  decreases,  r  must  increase.  The  following 
parameters  were  given  in  order  to  record  the  variation  of  r 
with  respect  to  H.  R=2000  km,  u=. 25,  E  =  7x1012,  f  =  1 0 0  bars, 
$=0.12  E-8  cm,  (table  2). 

This  table  clearly  shows  the  role  that  H  plays  in  the 
range  at  which  effects  due  to  a  particular  stress  drop  can 
be  observed. 

The  double  Benioff  zone  which  seems  to  be  present  in 
figures  23b,  24a  and  27  could  be  justified  by  some 
observations  in  some  island  arcs  (Engdahl  and  Scholz,  1977 
and  Uyeda,  1977).  The  presence  of  this  kind  of  weakness 
pattern  would  indicate  that  in  this  subduction  zone  the 
downgoing  slab  has  an  angle  of  about  45°.  However,  the  lack 
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THICKNESS 
OF  THE 
ELASTIC 
PLATE 
(H  IN  KM) 


RANGE  OF  THE 
OBSERVATION 
(r  IN  KM) 


20 

1876 

30 

1770 

40 

1641 

50 

1490 

60 

1318 

70 

1  118 

80 

888 

90 

610 

100 

200 

Table  2....  Shows  the  effect  of  the  thickness  of  an  elast 
plate  on  the  distance  from  an  applied  force  at  which  a 
displacement  z  takes  place. 
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of  this  pattern  in  figure  24b  makes  the  last  remark  weak  in 
the  sense  that  earthquakes  in  any  of  the  Benioff  zones 
should  have  aftershocks  in  both  of  these  zones  due  to  the 
proximity  between  them. 

To  interpret  figure  25  two  assumptions  made  should  be 
recalled.  First,  the  shear  energy  due  to  the  presence  of  the 
source  stays  almost  entirely  in  the  elastic  part  of  the 
model  and  secondly,  that  the  stability  distribution  does  not 
necessarily  represent  a  large  change  in  the  state  of  stress 
of  particular  locations,  but  rather  a  stability  estimate 
with  respect  to  the  rest  of  the  model  at  a  given  time.  With 
this  in  mind  figure  25b  shows  that  after  the  earthquake  the 
stability  of  the  elastic  part  has  decreased  to  values 
smaller  than  those  of  the  inelastic  part  which  is  supposed 
to  be  in  the  same  stress  state  as  in  figure  25a.  That  is, 
the  stable  part  of  figure  25a  remains  the  same  but  the 
unstable  part  has  become  more  stable. 

The  long  range  stabilization  of  figure  25  and  the  very 
localized  one  shown  in  figure  26b  are  in  accordance  with  the 
observations  that  at  shallow  angles  of  subduction  the 
effects  of  earthquakes  are  felt  at  longer  distances  than 
those  at  slabs  with  bigger  angles  (Kanamori,  1977  and 
Kelleher,  1974).  This  means  that  the  distribution  of 
probability  function  changes  due  to  external  forces  is 
greatly  dependent  on  the  geometry  of  the  subduction  model 
and  on  the  magnitude  of  the  frictional  forces  or  coupling 
between  the  two  plates.  For  a  bigger  coupling  or  a  smaller 
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subduction  angle  there  is  a  larger  effect  range.  ^ 

The  common  feature  in  all  the  after-earthquake  figures 
is  the  presence  of  an  instability  ring  around  the 
hypocenter.  This  represents  the  immediate  aftershock  area 
due  to  the  large  event.  This  ring  is  somewhat  vague  in  some 
cases,  as  in  figure  28b  where  only  an  arc  above  the  focus 
can  be  seen.  However,  since  the  earthquakes  are  simulated  as 
stress  drops  in  a  particular  location,  this  stress 
relaxation  should  stabilize  its  immediate  neighbourhood  for 
large  earthquakes. 

The  concentration  of  seismic  activity  in  the 
surroundings  of  the  focus  of  a  recent  large  earthquake  is  in 
accordance  with  the  secondary  stage  of  the  doughnut  pattern 
(Mogi,  1969),  which  states  that  with  the  occurrence  of  the 
great  earthquake  the  observed  seismic  pattern  of  a  region 
disappears  and  the  activity  is  concentrated  near  the  focal 
region.  This  kind  of  nearby  aftershock  area  is  also 
suggested  by  other  seismic  activation  patterns  observed 
( Kei 1 i s-Borok  et  al.,  1980. a). 

The  deep  high  probability  areas  shown  in  figures  27  and 
28  are  due  to  the  nature  of  the  approach.  A  quasi-static 
approach  of  the  problem  has  been  taken  here.  First,  the 
model  is  taken  at  equilibrium  and  at  a  given  time  the  forces 
are  set  free  to  act  in  the  model,  from  that  the  final 
distributions  of  probability  presented  were  calculated.  This 
kind  of  analysis  shows  for  a  viscous  material  its 
probability  of  yield.  It  does  not  assume  a  particular  mode 
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of  yielding.  ^ 

Although  the  mechanism  of  deep  and  intermediate 
earthquakes  is  stricly  unknown  this  kind  of  analysis  will 
only  be  valid  if  a  rapid  form  of  yielding  of  the  material  is 
possible . 

Several  probable  causes  have  been  proposed  for 
intermediate  and  deep  earthquakes.  Isacks  and  Molnar  (1971), 
suggested  that  if  the  slab  cannot  penetrate  beyond  650  km, 
then  this  involves  an  upward  transmission  of  compressional 
stress  that  could  cause  deep  earthquakes.  Griggs  (1972),  has 
shown  that  this  kind  of  resistance  would  result  in  some  kind 
of  buckling  of  the  subducting  plate,  for  which  there  is 
little  evidence. 

The  most  widely  accepted  theory  is  that  due  to  Ringwood 
(1975)  which  suggests  that  intermediate  and  deep  events  are 
caused  by  phase  transformations  of  the  material  that  is 
being  subducted,  which  result  in  sudden  volume  changes  in 
the  slab  that  might  create  instabilities.  However,  it  is  not 
very  well  understood  how  phase  changes  could  occur  fast 
enough  to  produce  earthquakes. 

Spanos  (1977),  suggested  that  intermediate  and  deep 
earthquakes  could  be  due  to  a  shear  heating  instability  in 
the  fluid  layer  between  the  slab  and  the  mantle  outside  the 
upper  boundary  of  the  subducting  slab.  He  demonstrated  that 
if  the  earthquake  occurs  within  this  layer  it  will  reproduce 
the  observed  focal  mechanism  for  deep  earthquakes.  That  is, 
the  principal  stresses  are  developed  within  the  downdipping 
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slab  and  in  the  direction  of  motion  of  the  platedsacks  and 
Molnar,  1971).  Sharp  changes  in  material  properties  have  not 
been  included  in  this  calculation.  This  has  resulted  in 
smooth  probability  curves  such  as  those  shown  in  chapter  3. 
In  this  particular  case  the  high  probability  zones  on  the 
downgoing  slabs  of  the  large  models  appear  without  gaps.  It 
is  not  the  purpose  here  to  find  a  mechanism  that  accounts 
for  the  presence  of  deep  earthquakes,  but  rather  to  see  the 
effects  of  a  sudden  stress  release  on  the  distribution  of 
probability  of  the  models,  and  how  these  distributions 
change  with  the  geometry  of  the  model  and  the  magnitude  of 
the  applied  forces. 

However,  the  sudden  displacements  at  great  depths  in 
the  form  of  earthquakes  is  controlled  by  the  amount  of 
heating  required  for  melting  and  by  how  well  defined  is  the 
melting  point  (Spanos,  1977).  If  the  temperature  of  the  slab 
is  far  from  the  melting  point  the  fluid  layer  will  be  narrow 
while  if  the  melting  point  is  well  defined  the  portion  of 
the  slab  that  reaches  this  temperature  will  have 
significantly  less  viscosity  than  the  rest  of  the  slab.  In 
any  of  these  two  cases  the  material  will  yield  almost  its 
entire  deformation  suddenly,  causing  total  melt  and  a 
complete  stress  relief.  This  kind  of  mechanism  generalizes 
our  results  to  the  ductile  regime. 
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4.5  Conclusions. 

The  thickness  of  the  elastic  part  of  the  lithosphere 
influences  the  range  at  which  long  range  precursors  can  be 
observed.  This  is  because  the  shear  energy  tends  to  remain 
within  this  elastic  layer,  that  is  the  viscous  material  acts 
as  a  barrier  and  rather  than  dissipate  this  energy  it  just 
does  not  permit  its  propagation.  Rupture  size  could  be 
linked  to  how  localized  an  anomalous  high  probabilty  is. 
Longer  areas  may  be  broken  in  more  high  probability  valued 
areas . 

Probability  values  are  relative  within  a  model 
calculation,  therefore  the  actual  value  is  just  an  indicator 
of  where  failure  is  more  likely  to  occur  for  a  set  of  given 
conditions.  It  is  more  reasonable  to  focus  this  kind  of 
analysis  on  the  study  of  the  distribution  of  high  and  low 
valued  areas  rather  than  on  the  values  themselves.  High 
probability  values  do  not  mean  that  an  earthquake  will 
necessarily  take  place,  but  just  that  at  a  particular 
location  failure  is  more  likely  than  at  other  parts  of  the 
same  model. 

The  decrease  of  probability  in  the  upper  tens  of 
kilometers  observed  in  figures  25b  and  26b  indicates  that 
the  whole  model  becomes  more  stable  at  these  depths  after 
the  occurrence  of  a  large  earthquake  and  that  aftershocks  in 
the  high  valued  zones  of  this  figure  are  the  most  likely 
kind  of  activity  to  occur  in  this  case. 
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The  result  that  after  an  earthquake  has  been  modelled, 
zones  of  high  risk  are  created  in  the  continental  plate 
presents  the  possibility  that,  if  this  plate  is  fractured 
and  partial  melting  occurs  at  depth,  then  this  molten  magma 
can  more  easily  reach  the  surface. 

No  difference  can  be  seen  in  the  range  of  effect  of  the 
deep  earthquakes  modelled  in  figures  27  and  28.  This  is  so 
because  they  are  both  originated  outside  the  elastic  parts 
of  the  plates. 

The  assumption  that  the  lithosphere  can  be  modelled  as 
an  elastic  layer  over  an  inelastic  one  seems  to  be  a 
plausible  one  and  the  solution  to  Boussinesq's  problem  and 
the  thin  plate  approach  shown  in  the  discussion  of  results 
add  credibility  to  this  kind  of  model. 

I  believe  that  Barenblatt  is  correct  by  suggesting  that 
local  concentrations  of  stresses  are  responsible  for  the 
presence  of  seismic  patterns  in  some  areas,  and  that  this 
can  be  modelled  by  adopting  Newman  and  Kn-opoff’s  ideas  of 
crack  populations  or  by  just  considering  the  presence  of 
geological  discontinuities  in  the  lithosphere.  However,  the 
range  of  effect  that  these  seismic  patterns  might  have  is 
strongly  linked  to  the  regional  characteristics  of  the 
material  in  which  these  concentrations  take  place.  That  is, 
it  is  controlled  by  the  thickness  of  the  elastic  part  of  the 
lithosphere . 
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5.  Earthquake  risk  studies  for  the  problem  of  induced 
seismicity  on  water  reservoirs  applied  to  the  Itzantun  site. 

This  chapter  is  an  application  of  the  instability 
function  described  in  section  2.4  to  a  site  located  in  the 
mountain  belt  associated  with  the  Middle  America  trench 
south  of  the  Tehuantepec  ridge,  the  Cocos-Car ibbean 
subduction  zone.  It  will  provide  some  insight  into  the  role 
of  this  system  in  the  local  seismicity  observed  in  the 
Itzantun  area.  This  part  of  the  thesis  deals  with  the  use  of 
time  dependent  instability  curves  to  make  estimates  of  the 
time  at  which  a  particular  location  becomes  more  unstable. 
Here  the  changes  of  the  stress  components  are  calculated 
analytically  by  means  of  the  formulation  described  in 
appendix  B. 

The  Itzantun  site  is  located  approximately  100  km 
inland  of  the  boundary  between  the  North-Amer ican  and  the 
Cocos  plates.  It  is  located  a  few  tens  of  kilometers 
North-West  of  the  portion  of  the  Motagua  Fault  that  serves 
as  boundary  between  the  North-Amer ican  and  the  Caribbean 
plates.  This  fault  is  the  westward  continuation  of  the 
tectonic  feature  known  as  the  Cayman  Trough.  As  mentioned 
earlier,  in  chapter  1,  one  of  the  features  of  plate 
tectonics  is  the  relation  between  the  mountain  belts  and  the 
plate  limits. 

In  the  portion  of  the  Middle  America  Arc  near  the 
Itzantun  area  seismic  activity  occurs  often.  It  is 
interesting  to  see  how  the  stability  distribution  behaves  in 
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this  particular  case  since  the  area  is  located  in  the  Sierra 
Madre  del  Sur.  This  is  the  mountain  belt  associated  with  the 
Middle  America  Subduction  Zone. 

Consolidation  theory  and  concepts  of  rock  failure  can 
be  used  to  evaluate  the  probable  risk  of  induced  seismicity 
as  a  result  of  filling  of  reservoirs.  The  evaluation  of 
seismic  instability  distribution  due  to  the  water  load 
superimposed  by  the  filling  of  a  reservoir  is  treated  in 
this  chapter.  As  an  example  of  a  particular  application,  the 
procedure  that  will  be  described  is  applied  to  a  particular 
reservoir.  This  evaluation  has  the  potential  to  indicate  the 
safest  way  to  fill  a  reservoir,  and  depends  only  on  the 
geometry  of  the  load,  the  rate  of  filling  and  the  geological 
structures  in  the  area.  The  stability  functions  as  discussed 
in  chapter  2  can  actually  be  a  measure  of  the  risk  of  having 
failure,  with  time,  for  a  particular  loading  history  in 
respect  to  a  plane  of  weakness. 

5.1  Evaluation  of  the  risk  of  induced  seismicity. 

Drawdowns  can  increase  the  risk  of  triggering 
earthquakes  in  this  area,  which  is  prone  to  thrust  faulting. 
It  is  possible  to  estimate  the  stresses  after  a  period 
during  which  the  water  level  is  maintained  and  a  decrease  in 
stresses  with  the  depth  of  the  observation  point. 

The  estimates  of  the  probable  induced  seismicity  are 
limited  as  the  residual  stress  in  the  area  prior  to  the 
impounding  is  unknown.  With  a  measure  of  the  residual 
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tectonic  stress  it  will  be  possible  to  determine  an  optimal 
filling  rate  to  reduce  the  probability  of  induced 
seismicity . 

The  relationship  between  mountain  belts,  continental 
limits  and  earthquakes  has  been  recognized  for  more  than  40 
years.  Gunnn  (1947)  suggests  that  the  parallel  distributions 
of  mountains,  ocean  deeps  and  volcanoes  follow  naturally 
from  the  properties  of  the  lithosphere,  if  it  is  assumed  to 
be  strong,  elastic  and  supported  by  an  underlying  weak 
asthenosphere .  He  proposed  that  the  generation  of  ocean 
deeps  were  possible  if  the  following  conditions  occur 
( figure  29 ) : 

1.  Down-folding  produced  by  horizontal  compression 
accompanied  by  a  thickness  of  the  sedimentary  layers. 

2.  Down-folding  derived  from  downward  convective 
circulation  of  the  underlying  magma. 

3.  Overthrusting  and  Underthrusting  at  shear  faults  in  the 
lithosphere. 

Gunn  (1947)  gives  some  ideas  that  have  become 
fundamental  for  plate  tectonics.  He  shows  that  in  order  to 
describe  the  gravity  anomalies  near  mountains,  deeps  and 
continental  boundaries,  the  lithosphere  cannot  be  weak  even 
if  moderate  stresses  are  applied  to  it  for  long  time.  He 
suggests  that  the  continental  lithosphere  must  be  at  a 
higher  state  of  stress  and  at  a  different  flexural  rigidity 
than  the  oceanic  one.  He  notes  that  most  earthquakes  in 
these  boundaries  represent  shear  fractures  and  that  their 
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Figure  29....  Calculated  deformation  of  a  shear  fractured 
lithosphere  subjected  to  a  horizontal  compressive  stress  of 
10*  dynes/cm2  (modified  from  Gunn,  1947). 


155 


spatial  distribution  in  depth  suggests  a  great  shear 
fracture  extending  throughout  the  surface  layers. 

5.2  General  concepts  on  induced  seismicity 

During  the  last  twenty  years,  it  has  been  observed  that 
large  engineering  projects  may  change  the  characteristics  of 
the  seismic  events  in  the  surrounding  region.  These  changes 
are  induced  by  changes  in  stress  that  are  a  result  of  man’s 
activities.  Among  the  activities  and  events  that  cause 
induced  seismicity  are  fluid  injection,  fluid  extraction, 
mining,  underground  detonations,  flooding,  and  reservoir 
impoundment  (Packer  et  al.  1977).  Here  we  will  deal  only 
with  reservoir  impoundment. 

There  are  many  examples  of  where  the  filling  of 
reservoirs  has  changed  the  character i st ics  of  events  in  an 
area.  These  changes  range  from  the  inducement  of  large 
magnitude  events  to  changes  in  micro-earthquake  activity. 

The  filling  of  large  reservoirs,  however  has  not  always 
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resulted  in  induced  seismicity.  Attempts  to  relate  induced 
seismicity  to  size  or  depth  of  a  reservoir  have  had  little 
success.  The  changes  in  seismic  activity  do  not  follow  a 
simple  pattern  (Gough,  1978).  Excellent  reviews  of  the 
observed  changes  have  been  prepared  by  Simpson  (1976),  and 
Gupta  and  Rastogi  (1976). 

Induced  seismicity  is  difficult  to  prove.  An  increase 
in  seismic  activity  in  areas  that  were  already  active  is 
difficult  to  attribute  entirely  to  the  filling  of  the 
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reservoir.  In  other  areas  the  pattern  of  seismic  events 
changes  radically,  and  there  seems  to  be  an  obvious 
association  with  the  filling  of  the  reservoir.  In  some 
areas,  there  appears  to  be  an  increase  of  seismic  activity 
during  the  initial  filling,  whereas  in  others,  the  increase 
occurs  some  years  after  filling.  There  appears  to  be  a 
correlation  between  the  water  depth  and  the  number  of 
earthquakes  at  some  reservoirs  (Withers  and  Nyland,  1978). 
And  there  also  appears  to  be  a  relation  to  the  rate  of 
filling  (Simpson  and  Negmatullaev ,  1981). 

The  amount  of  data  on  reservoir  induced  seismicity  is 
limited.  Up  to  1977,  there  had  been  55  reported  cases  of 
reservoir  induced  seismicity  (Packer  et  al,  1977).  Of  these, 
Packer  classifies  16  as  clear  cases,  35  as  questionable  and 
four  as  probably  not  reservoir  related.  Packer  et  al.  reach 
the  following  conclusions  regarding  induced  seismicity  due 
to  reservoir  loading: 

1.  The  initial  state  of  stress  in  the  ground  is  of  prime 
importance . 

2.  Failure  of  unfractured  material  as  a  result  of  reservoir 
filling  is  unlikely,  but  failure  is  likely  to  occur 
along  pre-existing  faults  in  fractured  material. 

3.  "Instantaneous"  stresses  generated  by  rapid  reservoir 
filling  lead  to  shear  stress  along  faults  without 
increasing  the  effective  stress. 

4.  Instability  along  faults  could  occur  at  great  depths  as 
shown  by  the  curvature  of  the  failure  envelope.  The 
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shearing  resistance  of  the  material  is  reduced  as  the 
confining  pressures  increase. 

There  is  by  no  means  unanimous  agreement  about  the 
existence  of  reservoir  induced  seismicity.  Other  authorities 
claim  only  3  clear  cases  of  reservoir  induced  seismicity 
exist.  The  difficulty  is  to  provide  a  viable  mechanism  for 
failure  caused  by  reservoirs  and  to  use  a  convincing 
stress-strain  relation  for  crustal  rock.  The  evidence  from 
other  reservoirs  indicates  convincingly  that  failure  can  be 
caused  by  relatively  small  external  influences  (i.e.  Gough, 
1978;  Wetmiller,  1981;  Raleigh  et  al.,  1976;  Pomeroy  et  al., 
1976;  Cook,  1976).  The  fact  that  statistically  rigorous 
observations  do  not  exist  for  reservoirs  does  not  imply  that 
reservoirs  can  not  induce  seismicity,  it  merely  means  that 
seismic  evidence  by  itself,  from  reservoirs  alone,  is  not 
sufficient  to  resolve  the  matter. 

A  consideration,  however,  of  the  existence  of  faults  on 
which  seismicity  is  known  to  occur,  the  fact  that  stimulated 
seismicity  has  been  observed  for  other  kinds  of  processes, 
and  the  fact  that  a  reasonable  physical  mechanism  for 
reservoir  induced  seismicity  can  be  postulated,  justifies 
modelling  studies  of  this  problem  to  determine  the  range  of 
risk. 

Any  prediction  of  seismicity  involves  assumptions  about 
the  stress-strain  relations  of  crustal  rock  and  the 
conditions  under  which  faults  will  fail.  The  largest  stress 
increment  due  to  large  reservoirs  is  of  the  order  of  10 
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bars.  Under  incremental  loads  of  10  bars  most  crustal  rocks 
deform  elastically.  Of  course  the  incremental  response  of  a 
rock  confined  under  103  bars  at  10  km  depth  may  be  different 
from  that  of  a  rock  at  the  surface,  but  its  elastic  nature 
remains  due  to  the  small  size  of  the  stress  increment. 
Therefore,  the  assumption  of  elastic  behaviour  is  plausible. 

The  variation  of  elastic  behaviour  can  be  deduced  from 
seismic  data.  Young's  modulus  deduced  from  seismic  data  for 
depths  from  0  to  25  km  varies  from  6x105  to  8x10s  kg/cm2. 
This  variation  is  small  compared  with  its  magnitude.  Hence, 
the  assumption  that  the  elastic  properties  are  constant  is 
reasonable  albeit  not  entirely  satisfactory.  Other 
authorities  (i.e.  Kirby,  1977;  Turcotte,  1974)  have 
considered  the  upper  25  kilometers  of  the  lithosphere  to  be 
elastic . 

Obviously  water  pressure  plays  a  crucial  role  in  the 
dynamics  near  a  reservoir.  The  simplest  extension  of 
elasticity  theory  that  takes  into  account  the  presence  of 
water  is  the  Biot  consolidation  theory.  It  is  normally 
applied  to  soils  and  is  justified  here  only  by  the  fact  that 
it  is  a  simple  tractable  extension  which  can  deal  with  the 
presence  of  pore  fluids  in  a  plausible  way.  It  may  not  be 
correct,  but  at  these  relatively  low  pressures  it  is  a 
reasonable  first  approximation. 

The  general  conclusion  from  the  observation  of  induced 
seismicity  is  that  reservoir  volume  is  not  always  a  reliable 
indicator  of  the  risk  of  induced  seismicity.  The  larger  the 
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volume,  the  greater  the  probable  risk,  but  there  is  always 
the  potential  for  suprises  such  as  was  encountered  at 
Hydro-Quebec  in  Canada  (Leblanc,  1978).  Manicouagan  3  on  the 
Canadian  Shield  caused  seismicity  changes  while  the  nearby 
Manicouagan  5,  twice  as  deep  and  with  a  considerably  larger 
volume,  has  not  induced  any  seismicity.  Manicouagan  3  has  a 
height  of  108  meters,  its  volume  is  1.04  x  1010  m3 . 

In  only  a  few  cases  have  the  depths  of  these  seismic 
events  been  determined  accurately.  Local  observations  and 
the  teleseismic  data  all  indicate  that  the  hypocentres  are 
shallow.  Gupta  et  al.  (1972)  have  determined  the  depths  and 
positions  from  the  events  at  Koyna  from  a  local  array,  and 
found  that  the  majority  of  the  events  occurred  at  a  depth  of 
less  than  10  km,  but  some  occurred  as  deep  as  30  km. 

Migration  of  seismic  events  has  also  been  observed  in 
some  reservoirs.  Simpson  (1976),  Simpson  and  Negmatullaev 
(1981)  and  Soboleva  and  Mamadaliev  (1976)  indicate  that  the 
events  at  Nurek  are  migrating  toward  the  reservoir. 

The  focal  mechanisms  (Bufe  et  al.,  1976.  Gough  and 
Gough,  1976  and  others),  observed  at  different  reservoirs 
are  consistent  with  the  types  of  pre-existing  faults  in  the 
neighbourhood.  At  Kariba,  Kremasta  and  Oroville,  dip-slip 
faulting  was  observed,  while  at  Koyna,  Hsinf engkiang  and 
Hoover,  the  mechanism  was  strike-slip  faulting.  At  Nurek  the 
induced  seismicity  is  occurring  along  a  series  of  thrust 
faults  connected  by  short  segments  that  show  strike  slip 
motion,  (Simpson  and  Negmatullaev,  1981).  Simpson  (1976), 
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Bell  and  Nur  (1978)  and  Withers  and  Nyland  (1976),  suggest 
that  rapid  lowering  and  raising  of  the  water  level  may  be  an 
important  factor  in  inducing  seismicity  in  regions  of  thrust 
faulting . 

The  magnitudes  of  the  main  shocks  near  reservoirs  have 
been  as  high  as  6.5  at  Koyna  (Gupta  et  al.,  1972),  6.3  at 
Kremasta  (Comninakis  et  al,  1968),  6.1  at  Hsenf engkiang 
(Wang  et  al.,  1976).  It  is  not  possible  to  give  an  upper 
limit  for  the  magnitude  of  induced  earthquakes,  as  the 
filling  of  reservoirs  acts  only  as  a  trigger  of  the 
pre-existing  stress. 

5.3  The  Itzantun  site 

The  Itzantun  site  is  in  the  state  of  Chiapas  in  the 
southern  part  of  Mexico  120  km  NE  of  the  city  of  Tuxtla 
Gutierrez  (figure  30)  It  is  in  a  region  with  several  rivers, 
the  most  important  of  which  are  the  Tlacotalpa,  the  San 
Pedro  and  the  Huitupan.  The  Tlacotalpa  flows  in  the  Itzantun 
gorge,  and  at  this  location  the  flow  is  2x10*  m3  of  water 
per  year.  The  geologic  formations  in  the  area  are  chiefly 
thick  assemblages  of  mudstones  and  massive  limestones. 

The  foundation  of  the  Itzantun  dam  will  be  sandstone, 
mudstone,  and  limestone  which  appears  reasonably 
homogeneous,  at  least  at  the  surface.  Many  fractures  in  the 
formations  near  the  dam  have  been  filled  with  calcite  but 
some  are  open  and  show  evidence  of  recent  movement  of  the 
order  of  centimeters. 
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Figure  30 
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Location  map  of  the  study  area 
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The  Itzantun  fault  crosses  the  reservoir  just  upstream 
from  the  dam  and  is  clearly  the  significant  structural 
feature  in  the  analysis  of  the  risk  of  induced  seismicity 
(figure  31).  The  gorge  itself  was  not  developed  along  a 
fault  zone.  The  irregular  directional  changes  of  the  river, 
and  the  fact  that  no  fault  breccia  was  found  in  a  borehole 
slanted  to  go  under  the  gorge,  indicate  that  the  river  has 
eroded  along  minor  fractures  and  joints.  Nevertheless  the 
strike  of  the  river  is  along  a  potential  failure  plane  (the 
known  faults  in  the  area  are  approximately  at  right  angles 
to  the  strike  of  the  river). 

5.4  Natural  seismicity  levels 

Historically  the  region  in  which  Itzantun  is  located 
has  been  the  origin  of  several  earthquakes.  Molnar  and  Sykes 
(1S69)  studied  21  earthquakes  in  the  portion  of  the  Middle 
America  Arc  that  lies  north-west  of  the  intersection  of  the 
arc  with  the  Motagua  fault,  and  in  general  found  that  the 
depth  of  the  events  increases  as  their  epicenters  move 
inland.  The  maximum  depth  determined  is  of  about  133  km; 
most  of  which  had  thrust  mechanisms.  Of  these  earthquakes 
the  majority  are  located  within  a  250  km  radius  of  Itzantun. 
Although  no  events  are  reported  by  Molnar  and  Sykes  in  the 
continental  part  of  the  North  Amer ican-Cocos  plate  boundary 
the  1976  Motagua  earthquake  of  magnitude  7.5  had  its  focus 
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Major  geologic  features  of  the  area  of 
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These  tectonic  characteristics  made  it  necessary  that 
the  project  of  building  a  deep  hydroelec t r ic  reservoir 
needed  an  assesment  of  the  natural  seismic  level  of  the 
area.  In  order  to  accomplish  this  a  nine  months' 
microearthquake  survey  was  carried  out  in  1979  in  the  area 
by  the  Institute  of  Engineering  of  The  National  University 
of  Mexico.  Six  portable  seismomenters  MQ-800  of  record  in 
smoked  paper  were  set  in  the  area  (figure  32).  More  than  70 
shallow  events  with  a  depth  <  30  km  were  detected  in  the 
area,  all  of  them  with  a  magnitude  <  2.2.  With  the 
information  obtained  in  this  first  period  of  recording  it  is 
not  possible  to  determine  if  any  of  the  faults  in  the  area 
is  active.  However,  the  majority  of  the  events  are  located 
in  the  north  of  the  intended  site,  which  is  the  part  where 
more  faults  are  present. 

The  local  network  described  was  able  to  record  during 
the  same  time  120  deep  earthquakes  with  a  magnitude  range 
from  3.0  to  5.0  (figure  33).  In  the  analysis  of  these  deep 
events  the  records  from  the  seismological  stations  Comitan 
(COM),  Chicoasen  (CSN)  and  Oaxaca  (VOH)  were  also  used.  The 
locations  of  these  earthquakes  (Novelo,  1978)  suggest,  that 
the  events  are  aligned  in  narrow  equal-depth  bands 
perpendicular  to  the  Tehuantepec  ridge  and  that  the  depth  of 
the  earthquakes  increases  with  the  distance  from  the  trench, 
which  is  in  accordance  with  Wilson's  (1963)  proposal  that 
the  direction  of  subduduction  may  be  determined  by  the 
aseismic  ridges.  The  hypocentral  locations  shown  in  figure 
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Figure  32....  Seismic  stations  and  epicentral 
determinations. 


166 


96°  95°  94°  93°  92°  91° 


A 

3.1 

4 

M 

< 

4.0 

LINE  •  DEPTH 

• 

4.1 

< 

M 

< 

5.0 

1- 115  KIV1 

2- 1  25 

3  -  135 
4-1  50 


Figure  33 ...  . 


Deep  events  located  by  the  Itzantun  network. 
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33  can  be  confined  to  a  certain  distance  from  the  trench 
(figure  34).  This  exercise  enables  one  to  fit  by  least 
squares  an  angle  of  subduction  of  45°  at  the  Cocos 
Cocos-North  American  plate  boundary. 


5.5  Model  studies 

As  a  first  approximation  it  is  possible  to  model  the 
problem  as  consolidation  of  a  water  logged  half-space.  Our 
computer  programs  for  2  and  3  dimensional  analysis  treat  the 
modelling  problem  by  considering  the  earth  to  be  a  uniform, 
isotropic  half-space  consisting  of  an  elastic  matrix 
affected  by  fluid  under  pressure.  This  material  is 
characterized  by  a  single  permeability,  a  relative  fluid 
matrix  compressibility,  a  coupling  factor  (or  hydraulic 
transmissibility )  for  the  bottom  of  the  reservoir,  and  two 
elastic  moduli.  The  reservoir  load  can  be  approximated  as  a 
"long"  2  dimensional  load  or  within  limits  treated  as  3 
dimensional. 

The  major  deficiency  of  this  approach  is  that  effects 
on  the  strength  of  faults  must  be  judged  qualitatively. 
Interpretive  examples  are  given  in  Withers  and  Nyland 
(1978).  They  point  out  that  the  incremental  stresses  due  to 
a  reservoir  are  rarely  large  enough  to  cause  failure  by 
themselves.  The  potential  for  failure  must  exist  and  may  be 
triggered  by  the  reservoir. 
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Figure  34....  Shows  the  deep  events  from  figure  33  in  a 
depth  vs  their  distance  to  the  Middle  America  trench 
diagram. 
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5.6  The  approach 

The  definition  of  stresses  in  porous  media  meets  with 
certain  difficulties,  but  some  heuristic  theory  has  been 
developed  to  deal  with  these  stresses.  Terzaghi ( 1 95 1 ) 
proposed  that  stresses  in  porous  media  are  a  'neutral 
stress',  the  stress  in  the  fluid  and  an  'effective  stress', 
the  difference  between  the  total  stress  prevailing  in  the 
fluid-filled  media  and  the  neutral  stress.  It  is  the 
effective  stress  that  causes  deformation  ( Scheidegger , 1 974 ) . 

Biot  (1941)  suggested  that  the  compaction  of  soils  is 
caused  by  a  phenomenom  called  'soil  consolidation'.  This 
means  that  the  settlement  is  caused  by  the  gradual 
adaptation  of  the  soil  to  a  load  variation.  Biot  made  the 
following  assumptions:  1)  Isotropy  of  the  material,  2) 
linear  stress-strain  relations,  3)  the  strains  in  the  media 
are  small,  4)  the  water  contained  in  the  pores  is 
incompressible  but  may  contain  air  bubbles,  and  5)  the  water 
flows  through  the  porous  skeleton  acording  to  Darcy's  law. 

With  these  assumptions,  Biot  developed  the  theory  for 
the  consolidation  of  porous  media,  the  basic  relations  that 
describe  the  phenomenom  are  given  by  Biot  in  a  series  of 
papers  published  since  1941.  Little  has  changed  in  this 
theory  since  then  (Rice  and  Cleary,  1976). 

In  order  to  approach  the  consolidation  problem  outlined 
earlier,  we  have  followed  the  technique  described  by  Withers 
and  Nyland  in  their  series  of  papers  (1976,  1977,  1978).  In 
order  to  solve  the  consolidation  equations,  use  is  made  of 
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the  displacement  functions  of  McNamee  and  Gibson ( 1 960 ) .  This 
implies  the  development  of  a  procedure  to  allow  us  to 
determine  the  double  Fourier  and  Laplace  transforms  of  the 
water  load.  The  Fourier  transforms  are  done  by  using  the 
Advanced  Mathematical  Library  of  the  array  processor 
(AP-190L),  which  allows  the  whole  computation  to  be 
overlapped  with  data  access  time.  This  permits  us  to  deal 
with  the  two  dimensional  transforms  of  the  load  at  a  given 
time  as  a  vector  and  determine  the  change  of  the  stresses  up 
to  that  time. 

Once  the  values  of  the  transforms  of  the  stresses  at 
the  desired  location  (corresponding  to  each  of  the  times  of 
the  known  load  history)  are  known,  there  is  adequate 
information  to  construct  a  curve  for  which  inverse  Laplace 
transform  will  give  the  behaviour  of  one  of  the  components 
of  stress  at  any  time  is  available.  From  these  components,  a 
failure  criterion  and  an  assumption  about  the  orientation  of 
a  plane  of  weakness,  we  calculate  estimates  of  stability  of 
a  point  in  the  formation.  The  inclusion  of  several  segments 
in  the  loading  history  curve  is  done  applying  the 
superposition  principle.  Thus  after  the  inverse  Fourier 
transform  is  performed  in  the  AP  for  a  given  component  of 
stress,  the  resulting  values  at  some  X,  Y,  Z,  are  the 
Laplace  transform  in  discrete  form  of  the  change  in  time  of 
one  component  of  stress.  The  result  is  a  function  of  time, 
which  defines  the  way  a  point  in  the  formation  moves  toward, 
or  away  from,  failure. 


■ 


i 


rr.i 

' 


171 


Instability  distributions  can  then  easily  be  determined 
by  using  the  procedure  described  in  section  2.4  of  this 
thesis.  The  primary  assumption  is  that  at  the  time  of 
impounding  the  reservoir  all  the  parts  of  the  porous 
halfspace  beneath  the  reservoir  were  stable.  That  is,  the 
instability  values  are  derived  from  the  incremental  stresses 
due  to  the  reservoir  and  the  geology  of  the  area,  and  they 
will  represent  how  instability  values  of  particular 
locations  are  changing. 

5.7  Discussion  of  Results 

I  attempt  here  to  evaluate  the  risk  of  induced 
seismicity  in  a  qualitative  way.  In  order  to  do  this  I  have 
made  and  justified  as  far  as  possible  a  number  of 
assumptions . 

1.  In  the  upper  25  km  of  the  earth  incremental  stress 
changes  cause  elastic  response  and  failure  occurs 
according  to  a  Mohr-Coulomb  failure  criterion. 

2.  The  in-situ  stresses  are  such  that  small  increments  on 
them  can  cause  failure. 

3.  The  effec.t  of  water  can  be  modelled  by  Biot 
Consolidation  theory. 

4.  A  uniform  halfspace  is  a  reasonable  approximation  to 
reality . 

5.  The  geologic  estimates  of  fault  orientation  define  the 
location  and  direction  of  expected  failure.  Intact  rock 
will  not  fail  under  reservoir  induced  loads. 
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Based  on  these  assumptions  the  results  shown  in  figures 
35  to  3^8  Vere  obtained.  The  stability  function  for  two 
loading  histories  consisting  of  monotone  increasing  loads, 
to  a  constant  load  were  calculated.  The  total  loading  times 
were  2  years  (figure  35)  and  16  years  (figure  36)  Both  have 
been  calculated  for  a  point  beneath  the  deepest  part  of  the 
reservoir  at  a  depth  of  1  km.  and  show  the  relation  of  the 
resulting  stress  to  the  rate  of  filling.  For  the  curve  where 
the  complete  load  is  reached  16  years  after  the  beginning  of 
the  impoundment,  the  increase  in  stress  is  much  smaller  than 
when  the  total  load  is  reached  after  2  years. 

Figure  35  also  illustrates  that  when  the  load  is  kept 
constant  for  a  certain  time  period,  the  stresses  begin  to 
decrease  to  a  limiting  value.  This  means  that  the  effect  of 
the  anomalous  stress  produces  changes  that  lead  to  an 
equilibrium  state  that  does  not  necessarily  have  to  be  the 
initial  state  of  stress  in  the  area.  This  can  be  thought  of 
as  related  to  the  existence  of  residual  stresses. 

The  risk  function  has  a  non-linear  dependence  on  the 
rate  of  filling  (table  3)..  This  table  shows  the  value  of  the 
stability  function  20  years  after  impounding  began  and 
maximum  value  attained  during  that  period.  This  is  done  for 
a  location  at  1  km  beneath  the  deepest  part  of  the  Itzantun 
reservoir.  The  loading  histories  used  to  obtain  this  table 
are  as  follows,  where  T  is  the  time  at  which  the  lake  was 
first  completely  filled,  and  D  is  the  maximum  depth  of  the 
reservoir.  After  .25  T  the  reservoir  had  water  up  to  .45  D. 
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Figure  35....  The  upper  graph  shows  a  loading  history,  where 
the  total  load  is  reached  after  2  years  of  loading.  The 
bottom  graph  shows  the  instability  at  1  km  depth. 
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Figure  36....  The  upper  graph  shows  a  loading  history,  where 
the  total  load  is  reached  after  16  years  of  loading.  The 
bottom  graph  shows  the  instability  at  1  km  depth. 


DURATION  OF 

STABILITY  FUNCTION  IN  BARS 

LOADING 

MAX. 

VALUE 

VALUE  AFTER 

IN  YEARS 

ATTAINED 

/ 

TIME 

20  YEARS 

2 

1 .54 

/ 

15 

yrs. 

1.21 

3 

0.95 

/ 

20 

yr s . 

0.95 

4 

0.51 

/ 

20 

yrs . 

0.51 

8 

0.15 

/ 

20 

yrs . 

0.15 

16 

0.04 

/ 

20 

yrs . 

0.04 

Table  3....  Variation  of  the  stability  function  for 
different  rates  of  filling  of  the  reservoir. 
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After  .5  T  it  had  .75  D.  At  .75  T  it  was  .9  D  full.  From  the 
time  T  the  reservoir  remains  filled.  It  is  obvious  that  a 
peak  in  the  stability  function  has  been  reached  during  this 
20  year  interval  only  for  the  first  case  in  table  1.  With  a 
faster  rate  of  filling  the  risk  of  reaching  failure  is 
higher,  for  the  first  case  the  risk  increases  sharply,  it 
reaches  its  maximum  value  2  years  from  the  beginning  of  the 
filling  of  the  reservoir  and  then  it  decreases  to  a  value  of 
about  .9  bars  and  remains  constant. 

The  rate  of  filling  of  the  reservoir  is  not  the  only 
way  in  which  artificial  lakes  could  change  the  seismic 
activity  of  an  area.  Some  changes  have  been  observed  after 
filling  and  draining  the  reservoir,  as  in  the  case  of 
Oroville  CA,  where  an  event  of  magnitude  5.9  occurred  after 
this  kind  of  loading  history  (Withers,  1977). 

In  order  to  see  the  effect  of  draining  of  a  reservoir 
on  the  instability  function  we  applied  this  kind  of  loading 
history  (figure  37).  This  example  shows  that  the  instability- 
function  for  unloading  tends  to  have  a  second  minimum,  in 
this  case  after  13  years.  The  effect  of  a  fast  decrease  in 
the  value  of  the  stability  function  must  generate  sudden 
changes  in  the  stresses  that  might  trigger  seismic  events. 

Other  histories  involving  reduction  of  loads  show  that 
if  lowering  the  water  level  is  done  rapidly,  the  negative 
slope  of  the  stability  function  moves  toward  a  horizontal 
positon.  This  reduces  the  risk  over  a  thrust  fault  but 
increases  it  for  a  normal  fault,  as  the  values  in  the  latter 
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Figure  37....  The  upper  graph  shows  a  loading  history  in 
which  an  unload  takes  place  before  the  total  load  is 
attained.  The  bottom  graph  shows  the  instability  at  1  km 
depth. 
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part  of  the  stability  curve  are  much  bigger. 

Decrease  in  loads  in  the  loading  history  result  in  a 
stability  decrease  that  attenuates  rapidly  with  the  depth  of 
the  observation  point.  For  a  depth  of  4  Km  the  effect  is  not 
observed  at  all  (figure  38). 

5.8  Conclusions  about  instability  at  Itzantun. 

Why  the  filling  of  some  reservoirs  causes  seismic 
events  is  poorly  understood.  We  give  no  firm  predictions  for 
Itzantun.  Studies  indicate  residual  stresses,  differences  in 
permeability,  and  differences  in  physical  properties  of  the 
formations  under  the  reservoir  may  determine  whether  there 
is  induced  seismicity  risk  or  not.  None  of  these  factors  is 
known  with  precision  at  Itzantun. 

The  changes  in  stability  around  a  water  reservoir  due 
to  the  presence  of  the  water  can  be  predicted  in  a 
qualitative  way  by  assuming:  a)  a  model  of  a  porous 
halfspace  consisting  of  an  elastic  matrix  saturated  with 
water,  b)  that  brittle  failure  can  occur  in  the  upper  ten 
kilometers  of  the  lithosphere  if  small  stress  changes  are 
made,  and  c)  that  the  effect  of  water  in  rocks  can  be 
approximated  with  Biot's  consolidation  theory. 

To  diminish  the  risk  of  induced  seismicity  at  Itzantun, 
the  filling  of  the  reservoir  should  be  as  slow  as  economics 
permits,  with  intervals  when  the  water  level  is  held 


constant . 
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Figure  38....  The  upper  graph  shows  the  same  loading  history 
as  that  of  figure  37.  The  bottom  graph  shows  the  instability 
at  4  km  depth. 
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It  is  unlikely  that  the  shallow  seismic  activity 
detected  in  the  Itzantun  area  could  be  somehow  artificially 
induced.  These  two  areas  are  also  different  because  of  the 
presence  of  deep  seismicity  near  Itzantun,  which  seems  to 
occur  at  depth  on  planes  parallel  to  the  Middle  America 
trench.  The  distribution  of  the  deep  seismic  observed  events 
along  the  trench  is  in  accordance  with  this  result. 

South  of  the  Tehuantepec  ridge  the  observed  relative 
velocity  of  the  plates  is  less  than  in  the  northern  part. 

The  angle  of  subduction  in  the  southern  branch  of  the  Middle 
America  trench  is  45°.  This  is  less  than  than  the  30°  angle 
deduced  for  the  northern  part.  Therefore,  the  deep 
seismicity  observed  in  the  Itzantun  area  must  be  related  to 
the  Cocos-Car ibbean  subduction  mechanism.  This  45°  angle 
implies  that  the  shallow  activity  observed  in  the  region  is 
not  related  with  deep  stress  drops  or  due  to  the  lower 
coupling  between  these  two  plates.  Probable  causes  for  it 
could  be  upgoing  magmas  and  more  likely  the  shallow  events 
could  be  related  to  the  presence  of  the  Motagua  fault  in  the 
vicinity  of  the  area. 


6.  Discussion  of  results  and  further  studies  proposed 

In  the  preceding  chapters  of  this  thesis  I  have 

developed  the  suggestion  of  treating  tectonic  problems  with 

stability  measures.  The  concepts  are  very  general  but  once 

they  are  understood  they  can  be  explored  by  making  changes 

in  the  parameters  involved.  These  changes  can  be  in  the 

definition  of  the  failure  envelope,  the  value  of  the  shear 

angle  used,  and  the  interdependence  of  the  material 

properties  of  the  models.  Changes  that  are  physically 

reasonable  do  not  usually  affect  the  nature  of  the 

instability.  In  some  cases  they  might  enhance  the 

distributions.  For  this  reason  it  is  important  to  justify  as 

much  as  possible  the  use  of  a  particular  set  of  parameters. 

The  choices  taken  throughout  this  thesis  are  not  necessarily 

optimal  ones  for  each  case,  they  are  in  some  cases 

simplistic  to  illustrate  the  method.  In  general  they  are 

representative  of  what  is  known  about  the  materials  involved 

in  each  model.  These  particular  choices  have  been  made  for 

* 

the  reasons  outlined  in  the  text  of  the  chapters  where  they 
were  applied. 

Following  the  same  procedure  many  other  studies  are 
possible.  One  of  the  most  interesting  ones  is  that  of  a 
triple  junction,  for  which  three  dimensional  analysis  is 
required.  Another  could  be  the  study  of  seismic  active  zones 
with  complicated  geology  in  order  to  refine  the  initial 
model  in  the  relation  to  the  observed  seismicity.  Geothermal 
zones  might  also  benefit  from  this  kind  of  treatment,  by 
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following  in  general,  the  procedure  outlined  in  chapter  IV. 

6.1  The  Cocos-North  Amer ica-Car ibbean  triple  junction. 

The  Middle  America  trench  is  the  arc  where  the  Cocos 
plate  subducts  beneath  the  North  American  and  the  Caribbean 
plates.  The  boundary  between  these  last  two  near  the  trench 
is  the  Motagua  fault,  which  is  the  continuation  of  the 
Cayman  trough.  These  three  plates  meet  in  an  unstable  triple 
junction  (figure  39).  The  area  surrounding  this  triple 
junction  is  the  object  of  this  study. 

Along  the  trench  several  small  seismic  gaps  exist  that 
correlate  with  the  topographic  highs  of  the  ocean  floor, 
which  are  associated  with  the  fracture  zones  and  the 
aseismic  ridges.  The  Tehuantepec  ridge  is  in  the  area  of 
interest  for  this  work.  It  is  probable  that  at  this  location 
slip  occurs  aseismically  or  that  seismic  slip  takes  place  in 
large  events  rare  in  time  or  perhaps  the  subduction  process 
has  been  severely  modified  (McNally  and  Minster , 1 98 1 ) . 

The  seismic  characteristics  of  the  area  have  been 
summarized  in  figure  39.  The  shallow  activity  (depth<  70  km, 
figure  39a),  occurs  in  bands  parallel  to  the  plate 
boundaries.  Taking  the  Middle  America  trench  as  divided  by 
the  triple  junction  (at  approximately  longitude  95°  W) ,  no 
significant  difference  can  be  observed  between  the 
North-West  and  the  South-East  parts  of  the  trench.  However, 
figure  39b  shows  that  the  distribution  of  deep  earthquakes 
(depth  >  140  km)  is  completely  one-sided.  Deep  seismic 
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Figure  39....  Tectonic  features 
recorded  in  the  triple  junction 
earthquakes  (depth  <  70  km)  and 
events  (modified  from  Molnar  and 


and  epicentral  locations 
area,  a)  Shows  the  shallow 
b)  the  intermediate  and  deep 

Sykes,  1969) 
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activity  is  concentrated  mainly  in  a  narrow  band  on  the 
continental  side  of  the  South-East  branch  of  the  trench. 

Deep  earthquakes  occur  where  the  Cocos  plate  subducts 
beneath  the  Caribbean  more  frequently  than  in  the 
Cocos-North  American  subduction  system. 

In  order  to  make  a  first  approach  to  this  problem  a  two 
dimensional  elastic  model  in  the  plan  of  the  area  was 
prepared  (figure  40).  In  this  simplistic  model  the  Cocos 
plate  pushes  the  North  American  and  the  Caribbean  plates  in 
normal  directions  to  the  boundaries  between  them  and  the 
Cocos  plate.  The  Caribbean  plate  is  taken  as  fixed.  The 
velocities  of  subduction  assumed  for  the  Cocos  plate  are  4.0 
cm/year  for  Mexico  and  2.3  cm/year  for  Central  America 
(after  McNally  and  Minster,  1981).  This  model  also  takes 
into  account  the  shear  stress  developed  between  the 
Caribbean  and  North  American  plates  along  the  Cayman  trough. 
Figure  40a  shows  this  model  in  which  all  the  plates  are 
taken  as  50  km  thick.  Due  to  the  nature  of  the  model  it  does 
not  permit  the  modelling  of  any  kind  of  three  dimensional 
structures  such  as  subduction  zones.  Although  it  is  highly 
idealized  the  earthquake  risk  distribution  (figure  40b) 
resulting  from  the  use  of  the  instability  function  as 
described  in  section  2.6  shows  the  following  interesting 
features . 

1.  It  predicts  a  banded  zone  of  high  risk  that  runs 

parallel  to  the  boundaries  between  the  Cocos  and  the 
other  two  plates  involved  in  the  system. 
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Figure  40....  Shows  the  two  dimensional  model  of  the  area  in 
the  upper  diagram,  and  the  resulting  instability  function 
for  this  model  in  the  bottom  one. 
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2.  This  risk  band  is  discontinuous  in  the  triple  junction 
area . 

3.  The  instability  varies  along  the  Cayman  trough  creating 
concentrated  areas  of  high  risk  along  it. 

It  is  remarkable  that  this  simple  model  delineates  some 
of  the  high  seismicity  areas  that  have  has  been  observed.  It 
also  predicts  the  non-uniform  seismic  behaviour  of  the 
Motagua  fault  in  which  events  have  been  rare  and  localized. 
However,  the  most  important  feature  that  this  model  predicts 
is  the  presence  of  the  Tehuantepec  gap,  which  is  a  feature 
that  has  been  the  object  of  many  investigations. 

The  agreement  between  the  results  of  this  model  with 
the  seismic  behavior  observed  for  the  region  suggests  the 
use  of  instability  functions  for  more  complicated  models. 
This  will  allow  a  better  understanding  of  the  tectonics  of 
the  area. 

The  three  dimensional  modelling  exercise  can  be 
approached  with  the  help  of  the  ADINA  program.  This  finite 
element  model  is  representative  of  the  geodynamical 
processes  that  are  thought  to  occur  in  the  area.  For  its 
understanding  the  discussion  of  the  results  presented 
throughout  this  thesis  is  of  primary  importance,  since  this 
model  involves  most  of  the  situations  that  have  been  treated 
in  the  different  parts  of  this  work.  The  characteristics  of 
the  model  shown  in  figure  41  are: 

1.  The  materials  of  the  upper  25  km  of  the  lithosphere  are 
treated  as  elastic. 
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2.  Visco-elasticity  is  assumed  beneath  this  depth.  The 
material  properties,  their  inter-dependence  and  their 
distribution  as  a  function  of  depth  is  taken  to  be  the 
same  as  in  the  models  of  the  different  slabs  on  chapter 
IV. 

3.  The  coupling  force  between  the  plates  varies  with  the 
angle  of  subduction  as  described  in  section  4.1. 

4.  At  the  North-Western  part  of  the  Middle  America  trench, 

the  Cocos  plate  subducts  beneath  the  North  America  plate 
at  an  angle  of  30°,  in  N37°E  direction  (McNally  and 
Mister ,  1981). 

5.  At  the  South-Eastern  branch  of  the  triple  junction  the 
Cocos  plate  moves  beneath  the  Caribbean  plate  at  a 
subduction  angle  of  45°  in  a  N26°E  direction  (McNally 
and  Mister ,  1981). 

6.  The  model  consists  of  550  nodes  and  405  parallelopiped 
elements  (figure  41). 

7.  Only  vertical  displacement  is  allowed  for  the  nodes 
located  in  the  vertical  outer  walls  of  the  model. 

8.  The  model  represents  an  area  of  2.75  x  10s  km2,  and  a 
depth  range  from  0  to  200  km. 

For  the  instability  study  of  this  model  the  probabilty 
function  described  in  section  2.6  was  used.  To  do  this, 
first  the  stress  components  for  each  of  the  nodes  of  the 
model  were  determined  in  the  same  way  as  in  chapter  IV.  The 
resulting  distribution  of  probability  of  seismic  failure  for 
this  model  is  shown  in  figures  42  and  43.  Horizontal  cross 


- 

N 

X- 


■ 


188 


12  345  6789  10 


Figure  41....  Shows  the  element  distribution  used  for  the 
three  dimensinal  model 
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sections  of  this  distribution  are  shown  at  a  depth  of  25  , 
50,  100,  and  150  km  (figure  42).  The  other  cross  section 
represents  the  probability  of  failure  along  the  plane  of  the 
boundary  between  the  plates  (figure  43)  that  is,  along  the 
upper  surface  of  the  subducting  slab. 

Figures  42  show  the  variation  of  the  distribution  of 
this  function  as  a  function  of  depth.  In  these  figures  the 
broken  dotted  line  between  points  A  and  B  represents  the 
upper  surface  of  the  subducting  slab.  The  other  dotted  line 
represents  the  vertical  projection  of  the  Cayman  trough  at 
depth . 

At  a  depth  of  25  km  (figure  42a)  the  high  probabilty 
area  is  a  wide  band  that  seems  to  follow  the  boundaries 
between  the  Cocos  and  the  other  two  plates.  The  interesting 
part  is  that  the  band  becomes  narrower  in  the  area  of  the 
Tehuantepec  gap.  This  high  probability  band  seems  to  reduce 
its  width  with  depth,  as  it  appears  as  a  narrow  band  in 
figure  42b  which  corresponds  to  the  50  km  cross  section. 
However,  it  remains  located  at  the  intraplate  boundary. 

Figure  42c  shows  that  at  a  depth  of  100  kilometers  the 
high  probability  area  has  vanished  completely.  This  figure 
also  shows  two  highly  localized  areas  where  the  probability 
function  attains  values  of  less  than  0.2,  one  of  them 
intersects  the  Cayman  trough  projection  and  the  other  is  by 
the  North-Western  part  of  the  Cocos-North  America  plate 
boundary . 
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Figure  42....  Shows  the  probability  function  distributions 
for  horizontal  cross-sections  at  depths  of  a)  25,  b)  50,  c) 
100  and  d)  150  km. 
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Figure  43....  Shows  the  probability  distribution  at  the 
upper  surface  of  the  slab  with  depth 
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This  low  probability  zone  seems  to  be  expanding  with 
depth  as  shown  in  figure  42d,  in  which  the  other  important 
feature  is  the  significantly  different  value  attained  by  the 
probability  function  at  both  sides  of  the  Cayman  trough  at 
depth.  An  area  of  probability  values  higher  than  .6  is 
present  inland  of  the  South-Eastern  branch  of  the  Middle 
America  trench.  All  these  can  be  interpreted  as  follows:  if 
deep  events  are  to  occur  they  should  be  focussed  to  the 
south-east  of  the  triple  junction.  This  is  in  accordance 
with  the  deep  events  described  in  chapter  V. 

The  other  point  to  be  made  from  figure  42  is  that 
although  the  probability  along  the  Cayman  trough  is  by  no 
means  uniform  with  depth,  little  can  be  said  that  relates  to 
the  localized  activity  observed  here.  It  seems  to  behave 
like  the  upper  part  of  the  North-Western  branch  of  the 
triple  junction.  Its  results  seem  anomalous  mostly  because 
both  areas  are  separated  by  the  Tehuantepec  gap  area. 

Figure  43  shows  that  the  distribution  of  the 
probability  function  is  discontinuous  across  the  zone 
located  just  North  of  the  triple  junction.  This  is  where  the 
Tehuantepec  gap  is  located.  It  confirms  the  increase  of 
stability  with  depths  inferred  from  figure  42,  and  also 
shows  a  big  difference  in  the  distribution  of  probability  at 
depth  between  both  sides  of  the  Caribbean-North  America 
boundary.  This  figure  indicates  little  about  the  behaviour 
at  shallow  depths  probably  because  of  the  size  of  the 
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6.2  Stress  determination  on  seismic  active  areas. 

In  regions  in  which  the  general  structure  of  the  upper 
crust  has  been  determined  to  some  extent,  and  in  which 
seismic  events  are  present,  instability  studies  can  be  used 
to  determine  the  state  of  stress  of  the  area  and  help  to 
refine  the  understanding  of  the  geologic  processes  which  are 
still  shaping  it. 

For  instance,  consider  a  thick  uniform  sedimentary 
layer  in  which  an  igneous  body  is  intruding  (figure  44). 

This  represents  the  instability  distribution  from  a  two 
dimensional  model  of  this  problem,  using  the  function 
described  in  section  2.3.  Figure  44a  shows  the  case  in  which 
the  intrusive  mechanism  is  considered  no  longer  active.  That 
is,  the  model  is  in  equilibrium.  Figure  44b,  on  the  other 
hand,  considers  that  the  intrussive  igneous  body  is  being 
pushed  toward  the  surface  by  a  vertical  force. 

From  the  comparison  of  both  distributions  of  figure  44, 
it  is  apparent  that  the  introduction  of  buoyancy  forces 
changes  the  distribution  of  risk,  in  such  a  way  that  the 
probability  for  the  occurrence  of  events  in  a  depth  range  of 
12  to  24  km  increases.  This  is  useful  in  two  ways.  One  is 
that  if  events  have  been  observed  at  particular  depths  the 
actual  buoyant  force  can  be  determined,  and  the  other  is  the 
opposite.  If  the  state  of  stress  can  be  determined  from 
other  sources  then  the  probable  depth  of  the  activity  of  the 
area  can  be  estimated. 
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Figure  44....  Shows  the  instability  of  a  model  of  a  magmatic 
intrusion,  a)  no  buoyancy  forces  are  applied.  b)Buoyancy 
forces  are  applied  at  the  bottom  of  the  intruding  body. 
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The  same  kind  of  approach  can  be  applied  to  problems 
involving  other  kinds  of  geological  structures  as  faults, 
folds,  horsts  or  grabens. 

6.3  Concluding  remarks. 

Instability  functions  are  useful  in  the  interpretation 
of  the  geodynamical  processes  that  take  place  in  the 
lithosphere.  The  results  are  highly  dependent  on  the 
accuracy  of  the  model  used  to  calculate  the  stresses.  This 
restricts  the  instability  distributions  to  be  representative 
only  in  a  regional  sense.  Local  high  resolution  studies  can 
be  done  if  a  detailed  model  is  available. 

The  results  obtained  here  for  the  Middle  America  region 
are  in  accordance  with  the  overall  knowledge  about  the 
tectonics  of  the  area.  The  observed  deep  seismicity  at 
Itzantun  and  the  probability  distributions  for  the  three 
dimensional  model  coincide  with  the  observed  seismicity. 
Shallow  events  (<  70  km  deep) ,  are  more  likely  to  occur  and 
deep  events  must  occur  more  often  south  of  the  Tehuantepec 
ridge.  On  the  Motagua  fault  the  probability  distribution 
indicates  that  events  should  occur  rarely. 

The  stress  drops  due  to  seismic  activity  to  the  north 
of  the  Tehuantepec  ridge  are  associated  with  a  shallower 
subduction  angle  than  those  occurring  in  the  southern  branch 
of  the  trench.  Therefore,  the  range  at  which  they  may  have 
any  effect  on  local  seismicity  of  other  areas  should  be 
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The  three  dimensional  probability  calculation  is  not 
able  to  predict  the  Tehuantepec  gap  in  a  satisfactory  way 
This  could  be  because  the  size  of  the  elements  used  to  model 
the  elastic  plates  is  too  big,  or  because  the  model  is 
inaccurate  for  this  part.  The  subducting  plate  is  modelled 
to  simulate  the  Cocos  plate  bending  from  a  angle  of  45°  in 
the  south  to  one  of  30°  in  the  northern  branch  of  the  trench 
across  the  Tehuantepec  ridge.  One  way  to  decrease  the 
probability  and  create  the  gap  is  to  assume  that  the  plate 
is  uncoupled  between  both  subduction  mechanisms,  which  is 
one  of  the  proposed  theories  for  the  existence  of  the  gap. 


. 
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Appendix  A.  Introduction  to  Finite  Elements. 


This  appendix  has  been  included  to  enable  the  reader  to 
follow  the  mathematical  procedure  followed  for  the 
determination  of  the  components  of  stress  in  chapters  3,  4, 
and  6.  It  shows  the  development  of  the  finite  element 
technique  for  its  most  simple  application  to  continuous 
elastic  materials  ( Zienkiewicz ,  1971).  The  extension  to 
other  kinds  of  materials  is  done  by  changing  the 
stress-strain  relationship  that  describes  the  behaviour  of 
the  material,  as  described  in  chapters  1  and  4  for 
visco-elastic  materials. 


Elastic  properties. 

For  a  three  dimensional  anisotropic  material  the  six 
stress  components  and  the  strain  components  in  a  Cartesian 
coordinate  system  are  related  by 
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principle  that  the  coeficients  of  the  above  matrix  must  be 
symmetrical  and  that  constants  are  sufficient  to  describe 
the  behaviour  of  the  material. 

In  geodynamical  applications  generally  and  in 
particular  in  rock  mechanics,  we  often  look  at  layered  or 
stratified  materials.  In  such  materials  symmetry  in 
behaviour  is  obtained  about  any  axis  perpendicular  to  the 
plane  of  stratification.  The  properties  in  this  plane  are 
independent  of  the  orientation  of  the  axis  in  the  plane. 

The  elastic  independent  constants  remaining  in  the 


stress-strain  relation  are: 
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(A. 2) 

By  analogy  with  the  isotropic  case  this  can  be  written 
as 

ex=(1/E1  )o  x  -  (v  i  /E  !  )  a  y  -  ( ,  /E  i  )  a  z 
e  y  =  (  v  i  /E  ,  )crx+(  \ /E2)o  y~(v  2/E2)o  z 
e  z=~(v  ,/E,)ox-(v2/E2)oY  +  ('\/E2)oz 

7  x  y  = (  1 /G 1  ) T x  y 
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7v  z  =  (  1/G2  )tv  z  =  {2(  1+y2  )/E2}rxy 

7  z  X  -  (  1 /G  t  )^xz 

(A. 3) 

if  we  limit  ourselves  to  the  solution  of  two  dimensional, 
plane  strain  problems.  Considering  plane  strains  parallel  to 
the  xy  plane  we  have 


e  z  ~ 7  y  z  — 7  x  z  —  0 


Substituting  this  into  the  above  we  get 


(A. 4) 


O  Z=V  !  (E2/E,  ) O  x  +  V 2  C7 y 


(A. 5) 


which  results  in  the  following  stress  strain  relations  in 
the  xy  plane: 


ex=(1/E1)(1-j'12E2/E1)ax-(j;1/E1)(1+i'2)ay 
e  y=“  (  V  ,/E  ,  )  (  1+y2)tfx+(  l/E2  )  (  1-*>22)c7y 
7  x  y  =  (  1 /G 1  )  T  x  y 

(A. 6) 

In  the  solution  by  stiffness  methods  it  is  more 
convenient  to  express  the  stresses  in  terms  of  strains. 
Writing  E2/Ei  =  V ,  we  have 

ax=E1(1+*'2)"1(1-J'2-277t'12)‘,{(1-t'22)ex+77t>1(1+i'2)ey} 

ay=E1(1+y2)‘,(1-t>2“277i'12)~,{i7J'1(1+j>2)ex+77(1-77*'12)ey} 

r  X  y =G 1 7  x  y 
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Or  in  matrix  form 

<?x  1  \ex  "| 

{ a  }  =  |  a  y  |  =  [  D  ]  |  e  y  |  =  [  D  ]  {  e  } 

/-vj  [7*VJ 

For  plane  stress  we  shall  have  by  taking 


(A. 7) 


(A. 8) 


C7Z— Ty2—  Tx2— 0 


the  following  stress  strain  relations 


(A. 9) 


ax=E  i  (  \-r\v  i  2  )  ‘  1  {  ex+^  ify) 
ay=E1  (  1-rjy  t  2)  '  1  {r?^1ex+r?ey} 

T x  y =G 1 7  x  y 

(A. 10) 

✓ 

The  Finite  Element  method. 

The  solution  of  any  stress  analysis  problem  can  be 
approached  via  the  formulation  of  the  governing  differential 
equations  or  by  the  use  of  various  energy  theorems. 
Essentially  the  finite  element  procedure  is  based  on  the 
latter,  and  uses  the  well  known  principle  of  minimum  total 
potential  energy.  If  U  denotes  the  strain  energy  expressed 
as  a  function  of  displacements  {6}  and  if  {P}  is  the  load 
system  associated  with  these  displacements  then  the  total 
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potential  energy  is 

0= U-U}  1  { P] 

(A. 1 1 ) 

and  has  to  be  minimized  in  the  region  considered.  Matrix 
notation  will  be  used  to  allow  for  any  load  system  applied. 

If  now  a  displacement  system  throughout  the  region  is 
specified  by  giving  displacement  values  at  each  of  the  nodes 
such  that  a  particular  node  affects  only  the  displacements 
in  adjacent  elements,  (figure  45),  then  an  approximate 
solution  can  be  obtained  by  differentiating  the  potential 
energy  with  respect  to  the  nodal  values,  element  by  element, 
and  adding  the  appropiate  differentials  of  all  elements  to 
minimize  the  total  potential  energy. 

In  any  element  the  displacements  {6}  at  any  point  can 
be  written  in  terms  of  the  nodal  displacements  { 5 d } 

(6}  =  [N]  { 6  d } 

(A. 12) 

from  which  the  strains  at  any  point  can  be  found  as 

{ e }  =  [B]  { 6  d } 

(A. 13) 

Then  by  Hooke's  law  the  stresses  are  given  as 

{  cr  }  =  [  D  ]  ({e}  —  {e0}) 


(A. 14) 
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Figure  45 ... . 
elements 


Subdivision  of  a  continium  into  finite 
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Now  the  strain  change  of  an  element  due  to  a  change 
A{5d}  may  be  expressed  as 


|"3Ud/96  j 


AUd  =  A{ 6  d } * |3Ud/96j |=JJA{e}t{a}d(vol) 


j^3Ud/36  k  J 


(A. 15) 


with  the  integration  taken  over  the  volume  of  the  element. 
Using  equations  (A. 13)  and  (A. 14) 

A(Ud)=A{6d}t{//[B]t[D][B]d(vol){6d}-/J[B]t[D]{e0}d(vol)} 

(A. 16) 

In  the  same  way,  if  the  external  loads  {P}  are  divided 
into  concentrated  loads  { C }  associated  with  nodal 
displacements  and  distributed  loads  {c},  the  variation  of 
the  loads  in  an  element  can  be  written  as 


A( - { 6 } T {P} ) =-A{ 6  d } t  | 3 { 6 } 'P/aS j 


={5d}t{-{C}-//[N]t{c}d(vol)} 


(A. 17) 


Thus  the  change  in  potential  energy  in  any  element  can 


be  expressed  in  terms  of  the  nodal  values  of  displacements 
associated  with  the  element  as 


’■ 
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| 90/3  6  j|  =  [K]{6d}_{C}-//[N]t{c}d(vol)-/i[B]t[D]{e0}d(vol) 


(A. 18) 


in  which 


[K]=//[B]t[D][B]d(vol) 


(A. 19) 


When  equating  the  individual  differentials  of  the  total 
potential  energy  to  zero  to  obtain  its  stationary  value,  the 
equations  corresponding  to  a  particular  node,  i,  are  set  up 
by  adding  all  the  appropriate  terms  K  f  j  6  j  and  the  ’i?  terms 
of  the  last  three  terms  of  equations  (A. 18)  for  each 
element . 

This  is  the  same  process  as  if  [K]  were  considered  as  a 
stiffness  matrix  of  a  structural  element, 

{C},{C’}=/J[N]t{c}d(vol)  and 

{C°}=//[B]t[D]{e0}d(vol)  (A. 20) 

as  the  loading  contributions  from  each  element  to  the 
equilibrium  of  a  node,  due  to  concentrated,  distributed  and 
initial  stress  forces  respectively.  Thus  the  resulting 
equation  at  a  point  ’ i '  due  to  equating  the  appropiate 
differentials  to  zero  will  then  be 


{Fj  ^ZCi+LC1  i+IC°  i  =  (IKi  , )61  +  (ZKi 2)62  +  - • • 


(A. 21 ) 
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The  summation  is  taken  throughout  all  the  elements. 

As  6j  affects  only  adjacents  elements  due  to  the  choice 
of  the  displacement  function  many  coefficients  in  the  right 
hand  side  will  be  zero  and  the  resulting  system  of  equations 

r«.i  ki 

[S] I  -  1  =  1 •  I 

MM 

(A. 22) 

will  have  a  narrow  band  form  which  will  be  found  of 
assistance  in  the  numerical  solution.  In  order  to  minimize 
this  band  width  it  is  important  to  number  the  nodes  across 
the  narrow  side  of  the  model  since  the  band  determines  the 
storage  space  (SS) 


SS=Nm 


(A. 23) 


were  m  is  the  maximum  nodal  difference  on  any  element  and  N 
is  the  total  number  of  nodes  in  the  model.  Therefore,  the 
number  of  operations  (NO)  that  have  to  be  carried  out  to 
determine  the  stiffness  of  the  model  is  of  the  order  of 


NO=Nm 2  (A. 24) 

Once  all  the  displacements  have  been  determined, 
stresses  in  any  element  can  be  calculated  by  equation 
(A. 14) . 
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Up  to  this  point  the  finite  element  description  has 
been  taken  for  a  general  case.  For  the  problem  of  a 
continuous  medium  divided  into  triangular  elements  (figure 
45),  the  displacement  throughout  the  element  with  components 
u  and  v  in  the  x  and  y  directions  can  be  defined  as 
{ 6 } =  fu] 


LI 

If  the  nodal  displacements  are  listed  as 

kl 


(A. 25) 


{  5  d  }  =  |  6  j 


L5kJ 


(A. 26) 


and  the  displacements  within  the  element  are  assumed  to  vary 


as 


kl 


{  6  }=  |~a  , 


+a2x+a3y | = | 1  x  y  0  0  0 


li- 
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(A. 27) 


then  the  six  constants  a  can  be  determined  by  equating  the 
six  nodal  displacement  components  with  x  and  y  coordinates 
taking  an  appropiate  value  as 
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| uk |  | 1  xk  yk  0  0  0  |  | a s  | 

00  0  1  xk  yk 


vk 


■JM 


by  equation  (A. 27) 

{5}=[l  x  y  0  0  o]  [A]"  ’  {6d}  =  [N]{6dJ 


L 


0  0  0  1  x  y 


(A. 28) 


(A. 29) 


The  strains  defined  as 

[ex  ]  [*9u/9x  ]  |"o 

{ e  }=  |  e y  |  | 9v/9y  |  =  |  0 

|^rxyJ  |^9u/9y+9v/9xJ  |^0 


1  0  0  0  o] 

0  0  0  0  1  |  [A] -  1 5d  =  [B] (6d } 
0  1  0  1  01 


(A. 30) 


define  the  [B]  matrix  of  equation  (A.  13). 

For  calculating  the  stresses  it  is  necessary  to  obtain 
the  D  matrix  in  the  x  and  y  coordinate  system.  In  the  form 
described  by  equations  (A. 9)  and  (A. 10)  the  matrix  is  valid 
for  a  system  in  which  the  x  is  perpendicular  to  the  strata. 
Calling  [ D T ]  the  matrix  in  which  a  load  system  of 


' 

- 
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coordinates  exists  defining 


(a’ }=[D' ]{e’ } 

\ 

where 

[DT  ]=E,  (  1+^2  )  '  1  (  1-^  2  -277^  ,  2  )  _  1  [(  \~v2 2  )  r\v  t(1+*>2)  0 


(A. 31 ) 


r]v  ,  (  }+v  2  )  7?  (  1-7^^  t  2  )  0  | 


0 


0 


G, 


J 


the  [D]  matrix  can  be  obtained  by  an  axis  transformation.  If 
the  angle  between  the  x'  and  the  x  axis  is  theta  as  shown  in 
figure  2, 

["cos2#  sin2#  -2sin#cos#~| 

{a}=|sin2#  cos2#  2sin#cos# | {a '  } 


sin#cos#  -sin#cos#  cos2#-sin2# 


Lsi 


=  [  T  ]  { a  '  } 

(A. 32) 

As  the  transformation  of  strains  follows  the  constant  work 
principle,  i.  e. 


£7  '  1  €  '  =0  t  e 

or  from  equation  (A. 32) 


(A. 33) 


{a}t([T]t)-1(e’}=at{e} 


. — ' 
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{e'  }  =  [T]  Me} 

{a' }  =  [ T ]  1 { a }  =  [ D T  ][T]t{e} 
{ a  }  =  [  T  ]  [  D  '  ]  [  T  ]  *  {  e  } 
therefore 


[D]=[T] [ D 1 ][T]1 

(A. 34) 

The  stiffness  matrix  given  by  equation  (A. 19),  is 
dependent  on  [D]  and  [B]  and  neither  of  these  matrices  is 
dependent  on  the  coordinates  of  integration,  thus  if  an 
element  of  unit  thickness  is  considered  in 


[K]=[B]< [D][B]A 


(A. 35) 

where  A  is  the  area  of  the  triangle  given  in  terms  of  the 
nodal  coordinates  as  the  determinant 


M  x  i  Y  .•  I 

A=1/2|1  xj  y  j  | 

|  1  xk  y k  | 

( A . 36 ) 

The  nodal  forces  and  body  forces  are  simply  computed 
from  equation  (A. 20)  with  all  the  matrices  there,  now 
defined.  For  the  element  described  the  body  forces  are 
assumed  constant  throughout  the  element  and  the  nodal  forces 
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are  simply  equal  to  total  body  force  acting  on  the  element 
assigned  in  three  equal  parts  to  each  of  the  nodes 
(Zienkiewicz  and  Cheung,  1965). 


. 


Appendix  B.  Calculation  of  stress  increments  in  porous 

media . 

The  description  of  stresses  in  porous  media  as 
mentioned  in  chapter  5,  deals  with  the  concept  of  effective 
stresses  (Terzaghi,  1951).  In  order  to  be  able  to  use  the 
instability  functions  it  is  necessary  to  calculate  the 
increments  of  these  stresses  when  a  change  in  the  water 
content  of  the  media  is  introduced.  This  involves  the 
assumption  that  before  any  variation  is  introduced  into  the 
medium  the  model  was  stable,  that  is,  no  seismic  movement 
occurs  until  the  variation  is  made.  The  effective  stress  is 
then  given  by  the  relation: 

a o  =  <7 n -P  (B  .  1  ) 

where 

a0=  effective  stress, 
an=  neutral  stress, 

P  =  pore  pressure. 

We  now  recall  the  assumptions  of  consolidation  theory 
(Biot,  1941),  which  were  described  in  chapter  5  of  this 
thesis , 

1.  Isotropy  of  the  material. 

2.  Linear  stress-strain  relations. 

3.  The  strains  in  the  media  are  small. 

4.  The  water  contained  in  the  pores  is  incompressible  but 
may  contain  air  bubbles. 
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5.  The  water  flows  through  the  porous  skeleton  according  to 
Darcy's  law.  (If  a  load  of  water  is  applied  over  a 
sponge  it  will  produce  a  gradual  settling,  depending  on 
the  rate  at  which  the  water  is  squeezed  out  of  the 
voids ) . 

Hooke's  equation  in  the  absence  of  fluid  can  be  written 

as : 

e  j  j  =  (  1  +v ) o j  j E" 1 -vo k k 6 j  j E" 1  (B.2) 

e  j  j  =  (9Ui/3Xj+9Uj/9Xj )/2 

where  v  is  Poisson's  modulus  and  E  is  Young's  modulus. 

If  P  is  the  pressure  of  the  water  contained  in  the 
pores  of  the  media,  and  H  describes  the  effect  of  the 
pressure  combined  with  the  previous  assumptions,  then  it  is 
possible  to  write. 

eij  =  (1+*')<7iji''1-j>crkk6ijE'1+P/3H6ii  (B.3) 

From  Biot  (1941)  we  can  write: 

8  = (ax x+ay Y+oz z ) /H i +P/R  (B.4) 

where  8  is  defined  as  the  increment  of  water  per  unit  volume 
of  the  material.  This  equation  following  Biot  (1941)  can  be 


written  as: 
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(7j  j  =2G  (  e  i  j  +  ( r?-  1  )  e  )  6  j  j-aa 

9  =  ari+o  /Q 

where : 

a=(  1  +i> )  G  {  3H  (  1-2*0  }"  1 

Q- 1 =R- 1 -a/H 

17  =  (  1  -v  )  /(  \-2v) 

G=E/ {2  +  2v ) 

e  =Vu 

a=(1-2^){2G(1-^)}-1 

* 

a i =a { 2G ( 1 }  1 

Introducing  this  notation  equation  (B.5)  gives 

G(V2Ui  +  (277-1  )ej  j  )-ao  i  j  =  0 

knowing  that  Darcy's  law  for  a  non  turbulent  flow  is: 


(B.5) 


(B.6) 


v i =-K9a/3x  j 


(B.7) 


\ 
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Where  K  is  the  permeability  of  the  material.  We  assumed 
that  the  water  is  incompressible,  therefore  the  increase  of 
water  content  in  an  element  of  the  solid  matrix  equals  the 
flow : 

90/3t=-3vx/3x-3vy/3y-3vz/3z  (B.8) 

With  (B.5),(B.6)  and  (B.8)  we  can  now  write 

KV2a-a3e/3t+Q" 1 3a/3t=0  (B.9) 

In  order  to  solve  equations  (B.6)  and(B.9)  we  will  use 
the  displacement  functions  of  McNamee  and  Gibson  (1960) 
following  the  method  used  by  Withers  (1976)  and  Withers  and 
Nyland ( 1976, 1978) . 

Taking  the  divergence  of  equation  (B.6) 

V2  (  2Grje-aa  )  =0 

Now  if  e  is  harmonic  (V2S=0)  it  follows 

2Grje-aa  =  2GdS/dz 

We  can  introduce  the  other  scalar  function  'E'  of 
Mcnamee  and  Gibson  (1960),  that  will  satisfy  this  equation 
if 

U=-3E/3x+z3S/3x 


V=-3E/3y+z3S/3y 


- 


■ 


\ 


\ 
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W=-9E/9z+z9S/9z-S 

P=2Ga '  1  (3S/3z-r?V2E) 

v 2  =-K2Ga "  1  (r?9e/9z  +  92S/3z2  ) 

e =V 2E  ( B  .  10) 

Substituting  (B.10)  into  equation  (B.9)  and  time 
scaling  t=tc,  we  get: 

V4E=d/dt (V2E)~2aid/dt (dS/dt )  (B. 1 1 ) 

and 

V2S=0  ( B . 12) 


where 

a  i  =G/ (  a  2Q+2Grj ) 

Equations  (B.11)  and  (B.12)  are  solved  by  Fourier 
transform  in  the  x,  y  variables  and  Laplace  transforming  in 
time  t.  The  notation  used  is 

Ei (pfq,z,t)=J /E(x,y,z,t)exp{i (px+qy ) } *dxdy  (B. 13) 

E2(p,q,z,s)=JE(p,q,z,t)exp(-st)*dt  (B.14) 


using 


■  . 
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d£2/dt=E 1(p,q,z,t=0)+s£2p,q,z,s)  ( B .  1  5 ) 

In  order  to  perform  the  Laplace  transformation  we  will 
now  introduce  the  following  procedure. 

The  depth  at  t  (ti<t<ti+1)  is 

z2(p,q,s)=/z1 (p,q, t ) exp (-st ) -dt 


=  z  2 ( p , q , t  j )  +  (t-tj ) {z2 (p,q,t i  +  ,) 


-z2 (p,q,t j ) } (t  j  +  i-t  s ) ‘ 1  ( B . 16) 

where  z2(p,q,t)  is  the  Fourier  transform  of  the  depth 
z(x,y,z).  The  Laplace  transform  in  t  of  this  function  is: 

z2 ( p ,  q ,  s ) =In - i z n - 1 {exp(-st n )~exp(-stn  _  !  )  } 


+(zn-Zn.i)(tn-tn-i)"{(tn- 1exp(“Stn . ! )“tnexp(-Stn ) )s‘ 1 

- (exp(-stn_exp(-st n - i ) )s" 1 }  +  s" 1 exp(-st „h)  (B  .  1 7 ) 
where:  s  is  the  Laplace  transform  parameter,  z,  is  the 
Fourier  transform  of  the  depth  at  the- time,  t„  is  the  nth 
time  at  which  a  depth  is  known,  and  m  the  total  number  of 
loading  points. 

We  now  have  a  formula  that  will  allow  us  to  determine 
z2(p,q,s)  once  the  Fourier  transforms  of  the  load  at  a 
finite  number  of  times  has  been  calculated.  These  Fourier 


transforms  are  done  by  using  the  advanced  mathematical 
library  subroutines  of  the  array  processor  (AP). 

The  AP-190L  is  a  high  speed  peripheral  arithmetic  array 
processor,  which  works  in  parallel  to  the  main  computer. 
Computation  of  fast  Fourier  tranforms  can  be  overlapped  with 
data  memory  access  time. 

And  as  indicated  in  chapter  5,  the  use  of  the  AP 
permits  us  to  deal  with  the  2  dimensional  transforms  of  the 
load  at  a  given  time  as  a  vector  and  determine  the  change  of 
stresses  up  to  that  time.  The  inverse  transforms  are  made 
without  using  complicated  algorithms  and  are  less 
time-consuming  in  the  computer  than  the  approach  followed  by 
Withers  (1977)  and  Nyland  And  Withers  (1976,  1978).  The  AP 
is  particularly  well  suited  to  the  large  numbers  of 
multiplications  and  additions  required  for  matrix  arithmetic 
due  to  its  highly  parallel  structure. 

Clearly  we  need  to  discretize  the  S  dependence  of  the 
problem.  We  evaluate  Laplace  transform  explicitly  at 
s  j  =  1 /t  j .  Thus  we  determine  the  parameters  A(p,q,s),  B(p,q,s) 
and  C(p,q,s)  involved  in  the  calculation  of  the  displacement 
functions  of  McNamee  and  Gibson  (1960).  From  these  the 
effective  stress  components  can  be  obtained. 

The  condition  imposed  by  the  use  of  the  Laplace 
transform,  that  we  start  calculations  at  time  t=0,  implies 
that  there  is  no  initial  compression  and  no  excess  pressure, 
this  means  that  the  result  will  reflect  the  effects  of  the 


anomalous  load  alone. 


. 
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For  t=0  we  have 


d/dz  Si ( t=0 ) =0  (B.18) 

Now  the  equations  (B.11)  and  (B.12)  simplify  to 

(-M+dVdz2 )  (-k2+d2/dz2  )£=-2aiSdS2/dz 


-K 2  +d  2S2/dz  2  =  0 


(B. 19) 


where  ju2=k2+S  and  k2=p2+q2 

The  solution  must  be  bounded  at  infinity  so  by  the 
method  of  variation  (Butkov,  1968,  p126)  we  may  write 

E= A  •  exp(  -kz  )  +B  •  exp(  -Atz  )  +ai  zC  •  exp(  -kz  ) 

S=C*exp-kz  (B.20) 

The  constants  A,  B  and  C  can  be  determined  from  the  boundary 
conditions:  at  z  =  0  a x  2  =a y  z  =  0 ,  and  at  the  surface  o z  z  is  the 
weight  of  the  water  alone,  which  is  represented  by  the 
compressive  load  -T2(p,q,0,s)  such  that: 

-T2/2G=kB(M-k)+Ck( 1-ai )  (B.21 ) 

Sustituting  E  and  S  into  equation  (B.10)  it  gives  the 
following  relations: 


• 


233 


aP/2G=-U2  +  k2  )B  +  2a  177-1  )  kC 

B=-(a«T2+aP2)  {2Gt?(m+X)  (ju+X)  }'  1 

=A{2 Gt?(m+X)  (m  +  X)  } '  1 

C=  (  2a  1 77-  1  )  -  1  k-  1  {aP2/2G+7?(M2-k2  )B} 

kA=~MB+a i C  ( B . 22 ) 

A=~ ( a  4T2  ~aP2 ) 

X=  ( 77+a «  )  kr? "  1 

3  n  —  ( 2a  1 77 —  1 )  ( 1 — a  1 ) 

The  excess  presure  in  the  pores  of  the  rock  is  then: 

P2/2G=Ck  (2ai77-1)exp(-kz)-7?B(ju2-k2)exp(-Mz)  (B.23) 

The  effective  stress  components  are  now: 

<7xx/2G=(92/9x2-V2  )E+z92S/9x2+9S/9z 

ayy/2G=(92/9y2-V2 )E+z9 2S/9y 2 +9S/3z 


oZ2/ 2G=( 32/9z2-V2 )£-z 9  2S/9z  2  +  9S/9z 


( B . 24 ) 


. 
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Once  the  stresses  have  been  determined  for  each  s  (  s 
is  the  Laplace  transform  variable)  by  means  of  the  procedure 
described  earlier  it  is  necessary  now  to  go  back  to  the  time 
domain,  the  inverse  Fourier  transforms  are  done  in  the  same 
way  as  the  direct  ones  but  in  order  to  perform  the  inverse 
Laplace  transform  we  used  the  following  algorithm. 

If  the  loading  history  can  always  be  represented  by  a 
curve  formed  by  several  segments  of  exponential  form,  then 
the  time  variation  in  any  of  the  components  of  stress  will 
also  be  represented  by  another  curve  consisting  of  segments 
of  exponential  form. 

Thus  once  we  have  the  values  of  the  transforms  of  the 
stresses  at  the  desired  location  for  each  s,  corresponding 
to  each  of  the  times  of  the  known  load  history.  We  actualy 
have  the  necessary  information  to  construct  a  curve  for 
which  inverse  Laplace  transform  will  give  us  the  behaviour 
of  one  of  the  components  of  stress  at  any  time. 

This  procedure  applied  to  the  three  principal 
directions  of  stress  will  allow  us  to  calculate  at  any 
moment  all  the  required  functions  to  develop  the  instability 
criteria  at  a  particular  location  with  respect  to  a  plane  of 
weakness . 

In  order  to  illustrate  let  us  consider  a  loading 
history  formed  by  only  one  segment.  The  curve  will  then  be 
represented  by 


f ( t ) =b{ 1 -exp( -at ) } 


( B . 25 ) 


. 

4  P 
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and  ’a'  and  ' b ’  are  determined  by  a  least  squares  technique 
from  known  values  of  f.  Then  the  transform  is 

f  (s)  =  /to{  l-exp(-at)  } exp (-st  )dt 

=b{s" 1 exp(-st ) + ( s+a) " 1 (exp(-st-at) }V 


F=f(s)=b{s“ 1  - ( a  +  s ) ~ 1 }  (B.26) 

When  (F) (s) ( s+a ) =b ( a+s~s ) =ba  we  can  write 
calling:  B=-ba 

B+aFs+s  2 F  =  0  ( B . 27 ) 

Given  some  values  of  F  it  is  routine  to  determine  B,  a 
and  hence  b. 

The  inclusion  of  several  segments  in  the  loading 
history  curve  can  be  done  applying  the  superposition 
principle  (depends  on  having  linear  constitutive  equations). 
Thus  after  the  inverse  Fourier  transform  is  performed  in  the 
AP  for  a  given  component  of  stress  the  resulting  values  at 
some  X,  Y,  Z,  are  the  Laplace  transform  in  discrete  form  of 
the  change  in  time  of  one  component  of  stress.  We  can  divide 
each  of  the  resulting  curves  into  exponential  segments  and 
determine  the  constants  a  and  b  that  represent  each  of  these 
segments.  Taking  the  inverse  Laplace  transform  of  each  of 
the  segments  involved  in  one  of  the  curves  and  adding  them 
up,  the  shape  of  the  variation  of  a  component  of  stress  can 
be  determined  for  all  times. 

The  proof  of  the  previous  statements  is  as  follows: 


:  - 
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If  f ( t ) i s  given  as : 

—  • 

f  ( t )  =b  * exp (-at )+c  t,<t<t2  (B.28) 

and  knowing  that  the  translation  of  a  Laplace  transform 
(Mathews  and  Walker,  1964)  by: 

/_{f(t-t1)}=exp(-t2s)L{f(t)}  (B.29) 

Then 

L( s)={exp(-t ! s ) -exp( -t 2 s ) }L{ f (t) } 

L(s)={exp(-t1s)-exp(-t2s)}(a+s)'1b  (B.30) 

Generalizing  (B.30)  for  t2  =  tn  and  t^t,,., 

L ( s )  =  {exp( -t  n - i s ) -exp( -t  n  s)}(an+s)',bn  (B.31) 

Therefore,  the  Laplace  transform  of  any  of  the  components  of 
stress  is  given  by  the  previous  equation  and  the  variation 
in  time  of  the  component  of  stress  is  given  as 

f n (t )=bnexp(-ant )+cn  tn-i^t<tn  (B.32) 

where  tn  is  the  nth  break  point  in  the  loading  history. 

To  test  the  previous  procedure  we  calculated  the  double 
Four ier-Laplace  transformation  of  several  shapes  of 
reservoirs  and  combined  them  with  different  loading 
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histories,  then  performed  the  inverse  transforms  in  the  way 
described  and  as  result  we  obtained  the  curve  that  we  fed  in 
as  the  loading  history.  For  example  the  resulting  curve 
(figure  46)  is  in  good  agreement  with  the  input  that 
originated  it  (Table  4). 


DEPTH  CIVII 
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TIME  CYR3 


Figure  46....  Resulting  curve  from  an  input  with  several 
exponential  segments 
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Depth  of  Water (mts . )  Time(yrs.) 

80  0.5 
110  1.0 
125  2.0 
140  2.5 
150  3.0 
155  4.0 
125  5.0 
140  5.5 
160  6.0 


Table  4....  Input  data  that  originated  Figure  46  using  the 
procedure  described  by  equations  (B.25)  to  (B.32). 
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